Ignore:
Timestamp:
05/28/04 11:56:06 (20 years ago)
Author:
rico
Message:
Changelog
Location:
trunk/MagicSoft/Mars/mtemp/mifae/programs
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc

    r4098 r4234  
    3636Float_t  lowdistcut=0;
    3737Float_t  uppdistcut=10;
    38 Float_t  sizecut=1800;
     38Float_t  lowsizecut=1800;
     39Float_t  uppsizecut=1800;
     40Int_t    nsizebins=1;
    3941UInt_t   onRate=1;
    4042Float_t  lowlgcut=0.01;
     
    130132  char evecut[256];
    131133  sprintf(evecut,"event%%%d==0",onRate);
    132   sprintf(gencut,"size>%f && dist>%f && dist<%f",sizecut,lowdistcut,uppdistcut);
     134  sprintf(gencut,"dist>%f && dist<%f",lowdistcut,uppdistcut);
    133135  cout << "General cut " << gencut << " (" << evecut << " for on-data)" << endl;
    134136
    135   // loop on widht and length cuts
    136   Float_t length=lowlgcut;
    137   while(length<=upplgcut)
    138     {
    139       Float_t width=lowwdcut;
    140       while(width<=uppwdcut)
    141         {
    142           TH1F* onhisto  = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90);
    143           TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90);   
    144 
    145           // define the width/length cut
    146           char varcut[256];
    147           sprintf(varcut,"length<%f && width<%f",length,width);
    148           cout << "Cutting " << varcut << endl;
    149          
    150           // define the on/off data cut
    151           char offcut[1024];
    152           sprintf(offcut,"%s && %s",gencut,varcut);
    153           char oncut[1024];
    154           sprintf(oncut,"%s && %s",evecut,offcut);
    155 
    156           // Fill the histos
    157           ton->Draw("alpha>>OnHist",oncut);
    158           toff->Draw("alpha>>OffHist",offcut);
    159          
    160           // Normalize histos
    161           const Int_t inibin = (Int_t)(frontiere/90.*nbins+1);
    162           Float_t level=0;
    163           Float_t leveloff=0;
    164           for(Int_t ibin = inibin; ibin<=nbins;ibin++)
     137  Bool_t isDiff = kFALSE;
     138  if(nsizebins<0)
     139    {
     140      nsizebins*=-1;
     141      isDiff=kTRUE;
     142    }
     143  else
     144    nsizebins+=1;
     145
     146  // loop on size, width and length cuts
     147  for(Int_t isb=0;isb<nsizebins;isb++) // loop on size
     148    {     
     149      cout << isb << endl;     
     150      Float_t minsize = isDiff ?
     151        TMath::Power(10,TMath::Log10(lowsizecut)+isb*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/nsizebins) :
     152        TMath::Power(10,TMath::Log10(lowsizecut)+isb*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/(nsizebins-1));
     153      Float_t maxsize = isDiff ? TMath::Power(10,TMath::Log10(lowsizecut)+(isb+1)*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/nsizebins) : 9999999999.;
     154     
     155      Float_t length=lowlgcut;
     156      while(length<=upplgcut) // loop on length
     157        {
     158          Float_t width=lowwdcut;
     159          while(width<=uppwdcut) // loop on width
    165160            {
    166               level+=onhisto->GetBinContent(ibin);
    167               leveloff+=offhisto->GetBinContent(ibin);
     161              TH1F* onhisto  = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90);
     162              TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90);   
     163             
     164              // define the width/length cut
     165              char varcut[256];
     166              sprintf(varcut,"size>%f && size<%f && length<%f && width<%f",minsize,maxsize,length,width);
     167              cout << "Cutting " << varcut << endl;
     168             
     169              // define the on/off data cut
     170              char offcut[1024];
     171              sprintf(offcut,"%s && %s",gencut,varcut);
     172              char oncut[1024];
     173              sprintf(oncut,"%s && %s",evecut,offcut);
     174             
     175              // Fill the histos
     176              ton->Draw("alpha>>OnHist",oncut);
     177              toff->Draw("alpha>>OffHist",offcut);
     178             
     179              // Normalize histos
     180              const Int_t inibin = (Int_t)(frontiere/90.*nbins+1);
     181              Float_t level=0;
     182              Float_t leveloff=0;
     183              for(Int_t ibin = inibin; ibin<=nbins;ibin++)
     184                {
     185                  level+=onhisto->GetBinContent(ibin);
     186                  leveloff+=offhisto->GetBinContent(ibin);
     187                }
     188              offhisto->Sumw2(); // needed to compute correct errors after normalization
     189              const Float_t norm = leveloff ? level/leveloff: 1;
     190              cout << "Normalizing by factor " << norm <<endl;
     191              offhisto->Scale(norm);   
     192             
     193              // compute significance/excess
     194              Float_t sig=0,bg=0,esig=0,ebg=0;
     195              Float_t significance=0;
     196              const Int_t signbins = inibin-1;
     197              for(Int_t ibin = 1; ibin<=signbins;ibin++)
     198                {
     199                  sig  += onhisto->GetBinContent(ibin);
     200                  esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
     201                  bg   += offhisto->GetBinContent(ibin);
     202                  ebg  += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
     203                }
     204              Float_t error= TMath::Sqrt(esig+ebg);
     205              significance = error ? (sig-bg)/error : 0;
     206             
     207              cout << "Excess: " << sig-bg << endl;
     208              cout << "N bkg: "  << bg << endl;
     209              cout << "Significance: " << significance << endl;   
     210             
     211              // save the result in file
     212              fout << minsize << '\t'
     213                   << maxsize << '\t'
     214                   << width << '\t'
     215                   << length << '\t'
     216                   << sig-bg << '\t'
     217                   << significance << '\t' << endl;
     218             
     219              delete onhisto;
     220              delete offhisto;
     221             
     222              width+=wdstep;
    168223            }
    169           offhisto->Sumw2(); // needed to compute correct errors after normalization
    170           const Float_t norm = leveloff ? level/leveloff: 1;
    171           cout << "Normalizing by factor " << norm <<endl;
    172           offhisto->Scale(norm);   
    173 
    174           // compute significance/excess
    175           Float_t sig=0,bg=0,esig=0,ebg=0;
    176           Float_t significance=0;
    177           const Int_t signbins = inibin-1;
    178           for(Int_t ibin = 1; ibin<=signbins;ibin++)
    179             {
    180               sig  += onhisto->GetBinContent(ibin);
    181               esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
    182               bg   += offhisto->GetBinContent(ibin);
    183               ebg  += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
    184             }
    185           Float_t error= TMath::Sqrt(esig+ebg);
    186           significance = error ? (sig-bg)/error : 0;
    187          
    188           cout << "Excess: " << sig-bg << endl;
    189           cout << "N bkg: "  << bg << endl;
    190           cout << "Significance: " << significance << endl;       
    191          
    192           // save the result in file
    193           fout << width << '\t'
    194                << length << '\t'
    195                << sig-bg << '\t'
    196                << significance << '\t' << endl;
    197          
    198           delete onhisto;
    199           delete offhisto;
    200          
    201           width+=wdstep;
    202         }
    203       length+=lgstep;
    204     }
    205  
     224          length+=lgstep;
     225        }
     226    }
    206227  fout.close();
    207228}
     
    260281        }
    261282
    262       // size cut
    263       if(strcmp(word.Data(),"SIZECUT")==0)
    264         ifun >> sizecut;
     283      // size cuts
     284      if(strcmp(word.Data(),"SIZECUTS")==0)
     285        {
     286          ifun >> lowsizecut;
     287          ifun >> uppsizecut;
     288          ifun >> nsizebins;
     289        }
    265290
    266291      // dist cuts
     
    317342  cout << "Off-data input file name(s): " << offFile << endl;
    318343  cout << "Output file name: " << outputFile << endl;
    319   cout << "Size cut: " << sizecut << endl;
    320344  cout << "Dist cuts (deg): [" << lowdistcut<< ","<<uppdistcut<< "]" << endl;
     345  cout << "Scanning size range: [" << lowsizecut << "," << uppsizecut << "]"<< " in " << nsizebins << " steps" << endl;
    321346  cout << "Scanning length range (deg): [" << lowlgcut<< ","<<upplgcut<< "], with step (deg): "<< lgstep << endl;
    322347  cout << "Scanning width range (deg): [" << lowwdcut<< ","<<uppwdcut<< "], with step (deg): "<< wdstep << endl;
  • trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizecuts.datacard

    r4117 r4234  
     1/////////////////////////////////////////////////////////
     2// Data card template for optimizeCuts executable      //
     3//                                                     //
     4// *** Compulsory cards are:                           //
     5// ONFILES, OFFFILES, OUTFILE                          //
     6//                                                     //
     7// the rest are optional                               //
     8// *** Do NEVER add a datacard with no value after it  //
     9/////////////////////////////////////////////////////////
    110
    211// Maximun number of on and off events to be processed)
     
    413
    514// On data acceptance rate (e.g 4 for taking 1/4 of whole data sample)
    6 ONRATE 2
     15ONRATE 3
    716
    817// Input file name pattern (On data)
    9 ONFILES  /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/crab20040215OnRotateCalA-D.root
     18ONFILES  /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OnRotateNoCalB.root
    1019
    1120// Input file name pattern (Off data)
    12 OFFFILES  /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/crab20040215OffRotateCalA-H.root
     21OFFFILES  /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OffRotateA-C.root
    1322
    1423// output file name
    15 OUTFILE  ./optimizeCrab_cal_fine.out
     24OUTFILE  ./prueba.out
    1625
    17 // Preliminar cuts (size in size units, distance in deg)
    18 SIZECUT  2000
    19 DISTCUTS 0.2 1.1
    2026
    21 // Length initial, final and step values
    22 LENGTHCUTS 0.15 0.35 0.005
     27// Preliminary distance lower and upper cuts (in deg)
     28DISTCUTS 0.3 1.2
    2329
    24 // Width initial, final and step values
    25 WIDTHCUTS 0.05 0.15 0.005
     30// Size initial, final and number of steps
     31// it will take logarithmic constant intervals
     32// positive number of steps means lower cut (integral)
     33// negative number of steps means lower and upper cuts (diferential)
     34SIZECUTS  500 10000 3
     35
     36// Length initial, final and step values (in deg)
     37LENGTHCUTS 0.15 0.35 0.1
     38
     39// Width initial, final and step values (in deg)
     40WIDTHCUTS 0.05 0.15 0.1
    2641
    2742
Note: See TracChangeset for help on using the changeset viewer.