loon: BinningScheme.cxx Source File

00001 #include "BinningScheme.h"
00002 
00003 #include "Configuration.h"
00004 
00005 #include "TSystem.h"
00006 #include "TDirectory.h"
00007 #include "TFile.h"
00008 #include "TH1.h"
00009 
00010 #include <iostream>
00011 #include <cmath>
00012 #include <cassert>
00013 
00014 using namespace OscFit;
00015 
00016 ClassImp(BinningScheme)
00017 
00018 static BinningScheme* fgBinningScheme = 0;
00019 
00020 BinningScheme* BinningScheme::Instance()
00021 {
00022   if( !fgBinningScheme ){
00023     fgBinningScheme = new BinningScheme();
00024   }
00025 
00026   return fgBinningScheme;
00027 }
00028 
00029 BinningScheme::BinningScheme()
00030 {
00031   fBinningScheme = kUnknownScheme;
00032 
00033   fBinningSchemeIsLoaded = 0;
00034 }
00035 
00036 BinningScheme::~BinningScheme()
00037 {
00038   std::map<Binning_t,Binning*>::iterator myIter = fBinningMap.begin();
00039 
00040   for( ; myIter!=fBinningMap.end(); myIter++ ){
00041     Binning* binningObject = myIter->second;
00042     if( binningObject ) delete binningObject;
00043   }
00044 }
00045 
00046 void BinningScheme::SetBinningScheme( BinningScheme_t theBinningScheme )
00047 {  
00048   if( fBinningSchemeIsLoaded ) {
00049     std::cout << "  warning: binning scheme is already loaded " << std::endl;
00050     return;
00051   }
00052 
00053   fBinningScheme = theBinningScheme;
00054 }
00055 
00056 Binning* BinningScheme::GetBinning( Binning_t binningType )
00057 {
00058   return Build( binningType );
00059 }
00060 
00061 Binning* BinningScheme::GetBinning( Run_t run, Containment_t containment, Flavour_t flavour, Resolution_t resolution )
00062 {
00063   return GetBinning( LookupBinning( run, containment, flavour, resolution ) );
00064 }
00065 
00066 Binning_t BinningScheme::LookupBinning( Run_t run, Containment_t containment, Flavour_t flavour, Resolution_t resolution )
00067 {
00068   switch( fBinningScheme ){
00069     case kScheme0:         return LookupBinningUsingScheme0( run, containment, flavour, resolution );
00070     case kScheme1:         return LookupBinningUsingScheme1( run, containment, flavour, resolution );
00071     case kScheme1A:        return LookupBinningUsingScheme1( run, containment, flavour, resolution );
00072     case kScheme2:         return LookupBinningUsingScheme2( run, containment, flavour, resolution );
00073     case kScheme2A:        return LookupBinningUsingScheme2( run, containment, flavour, resolution );
00074     case kCounterForMINOS: return kSingleBin;
00075     case kCounterForLBNE:  return kSingleBin;
00076     default:               return kUnknownBinning;
00077   }
00078 }
00079 
00080 Binning_t BinningScheme::LookupBinningUsingScheme0( Run_t run, Containment_t containment, Flavour_t flavour, Resolution_t resolution )
00081 {
00082   // Binning scheme for MINOS two-flavour analysis
00083   // =============================================
00084 
00085   // Atmospheric Neutrinos
00086   if( run==kAtmos ){
00087     if( resolution==kLoRes ){
00088       return kSingleBin;
00089     }
00090     else{
00091       return kAtmosLogLE;
00092     }
00093   }
00094 
00095   // Beam Neutrinos
00096   else{
00097     if( flavour==kNue ){
00098       return kBeamNUE;
00099     }
00100     else{
00101       if( containment==kRockMuon ){
00102         return kBeamRAF;
00103       }
00104       else{
00105         if (IsMINOSPlus(run)) return kBeamCCMinosPlus;
00106         return kBeamCC;
00107       }
00108     }
00109   }
00110 
00111   return kSingleBin;
00112 }
00113 
00114 Binning_t BinningScheme::LookupBinningUsingScheme1( Run_t run, Containment_t containment, Flavour_t flavour, Resolution_t resolution )
00115 {
00116   // Binning scheme for MINOS three-flavour analysis
00117   // ===============================================
00118  
00119   // Atmospheric Neutrinos
00120   if( run==kAtmos ){
00121     if( resolution==kLoRes ){
00122       if( Configuration::Instance()->DoingLowResolutionEnergyBins()==true )
00123         return kAtmosEnergy;
00124       else
00125         return kAtmosSingleBin;
00126     }
00127     else{
00128       if( GetContainmentType(containment)==kContainedVertex ){
00129         if( Configuration::Instance()->UsingAtmosLegacyBinning()==true ) 
00130           return kAtmos2D;
00131         else
00132           return kAtmos2DLoRes;
00133       }
00134       else{
00135         if( Configuration::Instance()->UsingAtmosLegacyBinning()==true ) 
00136           return kAtmosAngle;
00137         else
00138           return kAtmosAngleLoRes;
00139       }
00140     }
00141   }
00142 
00143   // Beam Neutrinos
00144   else{
00145     if( flavour==kNue ){
00146       return kBeamNUE;
00147     }
00148     else{
00149       if( containment==kRockMuon ){
00150         return kBeamRAF;
00151       }
00152       else{
00153         if (IsMINOSPlus(run)) return kBeamCCMinosPlus;
00154         return kBeamCC;
00155       }
00156     }
00157   }
00158 
00159   return kSingleBin;
00160 }
00161 
00162 Binning_t BinningScheme::LookupBinningUsingScheme2( Run_t run, Containment_t containment, Flavour_t, Resolution_t resolution )
00163 {
00164   // Binning scheme for LBNE three-flavour analysis
00165   // ==============================================
00166 
00167   // Atmospheric Neutrinos
00168   if( run==kAtmos ){
00169     if( resolution==kLoRes ){
00170       if( Configuration::Instance()->UsingAtmosLegacyBinning()==true ) 
00171         return kAtmosEnergy;
00172       else
00173         return kAtmosAngleHiRes;
00174     }
00175     else{
00176       if( GetContainmentType(containment)==kContainedVertex ){
00177         if( Configuration::Instance()->UsingAtmosLegacyBinning()==true ) 
00178           return kAtmos2D;
00179         else
00180           return kAtmos2DHiRes;
00181       }
00182       else{
00183         if( Configuration::Instance()->UsingAtmosLegacyBinning()==true ) 
00184           return kAtmosAngle;
00185         else
00186           return kAtmosAngleHiRes;
00187       }
00188     }
00189   }
00190 
00191   // Beam Neutrinos
00192   else{
00193     if (IsMINOSPlus(run)) return kBeamCCMinosPlus;
00194     if( Configuration::Instance()->UsingBeamLegacyBinning()==true ) 
00195       return kBeamCC;
00196     else
00197       return kBeamHiRes;
00198   }
00199 
00200   return kSingleBin;
00201 }
00202 
00203 Bool_t BinningScheme::CheckBinning( Run_t run, Containment_t containment, Flavour_t flavour, Resolution_t resolution, Int_t numBinsX, Int_t numBinsY )
00204 {
00205   Binning* myBinning = GetBinning( run, containment, 
00206                                    flavour, resolution );
00207 
00208   assert( myBinning );
00209 
00210   if( myBinning->GetNbinsX()==numBinsX 
00211    && myBinning->GetNbinsY()==numBinsY ){
00212     return true;
00213   }
00214   else if( myBinning->GetNbinsX()==1 ){
00215     return true;
00216   }
00217   else{
00218     std::cout << "  <warning> Input templates don't match binning scheme (" << AsString(fBinningScheme) << ")" << std::endl;
00219     std::cout << "            (run, containment, flavour, resolution) = (" << AsString(run) << "," << AsString(containment)
00220               << "," << AsString(flavour) << "," << AsString(resolution) << ")" << std::endl;
00221     std::cout<<myBinning->GetNbinsX()<<" "<<numBinsX<<std::endl;
00222     return false;
00223   }
00224 }
00225 
00226 Binning* BinningScheme::Build( Binning_t binningType )
00227 {
00228   Binning* binningObject = Lookup( binningType );
00229 
00230   if( !binningObject ){
00231     switch( binningType ){
00232       case kSingleBin:        binningObject = BuildSingleBin();         break;
00233       case kBeamCC:           binningObject = BuildBeamCC();            break;
00234       case kBeamCCMinosPlus:  binningObject = BuildBeamCCMinosPlus();   break;
00235       case kBeamNUE:          binningObject = BuildBeamNUE();           break;
00236       case kBeamRAF:          binningObject = BuildBeamRAF();           break;
00237       case kBeamHiRes:        binningObject = BuildBeamHiRes();         break;
00238       case kAtmosSingleBin:   binningObject = BuildAtmosSingleBin();    break;
00239       case kAtmosLogLE:       binningObject = BuildAtmosLogLE();        break;
00240       case kAtmosEnergy:      binningObject = BuildAtmosEnergy();       break;
00241       case kAtmosAngle:       binningObject = BuildAtmosAngle();        break;
00242       case kAtmosAngleLoRes:  binningObject = BuildAtmosAngleLoRes();   break;
00243       case kAtmosAngleHiRes:  binningObject = BuildAtmosAngleHiRes();   break;
00244       case kAtmos2D:          binningObject = BuildAtmos2D();           break;
00245       case kAtmos2DLoRes:     binningObject = BuildAtmos2DLoRes();      break;
00246       case kAtmos2DHiRes:     binningObject = BuildAtmos2DHiRes();      break;
00247       default:                binningObject = BuildSingleBin();         break;
00248     }
00249   }
00250 
00251   if( binningObject ) Load( binningType, binningObject );
00252   
00253   return binningObject;
00254 }
00255 
00256 void BinningScheme::TestBinningScheme()
00257 {
00258   Binning* binningObject = 0;
00259 
00260   binningObject = GetBinning( kSingleBin );
00261   if( binningObject ) binningObject->PrintBinningScheme();
00262 
00263   binningObject = GetBinning( kBeamCC );
00264   if( binningObject ) binningObject->PrintBinningScheme();
00265   
00266   binningObject = GetBinning( kBeamCCMinosPlus );
00267   if( binningObject ) binningObject->PrintBinningScheme();
00268   
00269   binningObject = GetBinning( kBeamNUE );
00270   if( binningObject ) binningObject->PrintBinningScheme();
00271 
00272   binningObject = GetBinning( kBeamRAF );
00273   if( binningObject ) binningObject->PrintBinningScheme();
00274 
00275   binningObject = GetBinning( kBeamHiRes );
00276   if( binningObject ) binningObject->PrintBinningScheme();
00277   
00278   binningObject = GetBinning( kAtmosSingleBin );
00279   if( binningObject ) binningObject->PrintBinningScheme();
00280 
00281   binningObject = GetBinning( kAtmosLogLE );
00282   if( binningObject ) binningObject->PrintBinningScheme();
00283   
00284   binningObject = GetBinning( kAtmosEnergy );
00285   if( binningObject ) binningObject->PrintBinningScheme();
00286 
00287   binningObject = GetBinning( kAtmosAngle );
00288   if( binningObject ) binningObject->PrintBinningScheme();
00289 
00290   binningObject = GetBinning( kAtmosAngleLoRes );
00291   if( binningObject ) binningObject->PrintBinningScheme();
00292 
00293   binningObject = GetBinning( kAtmosAngleHiRes );
00294   if( binningObject ) binningObject->PrintBinningScheme();
00295   
00296   binningObject = GetBinning( kAtmos2D );
00297   if( binningObject ) binningObject->PrintBinningScheme();
00298 
00299   binningObject = GetBinning( kAtmos2DLoRes );
00300   if( binningObject ) binningObject->PrintBinningScheme();
00301 
00302   binningObject = GetBinning( kAtmos2DHiRes );
00303   if( binningObject ) binningObject->PrintBinningScheme();
00304 
00305   return;
00306 }
00307 
00308 void BinningScheme::Load(Binning_t binningType, Binning* binningObject)
00309 {
00310   if( binningObject==0 ) return;
00311 
00312   std::map<Binning_t,Binning*>::iterator myEntry = fBinningMap.find(binningType);
00313   
00314   if( myEntry==fBinningMap.end() ){
00315     fBinningMap.insert( std::pair<Binning_t,Binning*>(binningType,binningObject) );
00316     fBinningSchemeIsLoaded = 1; // Binning Scheme is now loaded
00317   }
00318 
00319   return;
00320 }
00321 
00322 Binning* BinningScheme::Lookup(Binning_t binningType)
00323 {
00324   Binning* binningObject = 0;
00325 
00326   std::map<Binning_t,Binning*>::iterator myEntry = fBinningMap.find(binningType);
00327 
00328   if( myEntry != fBinningMap.end() ){
00329     binningObject = myEntry->second;
00330   }
00331 
00332   return binningObject;
00333 }
00334 
00335 
00336 Binning* BinningScheme::BuildSingleBin()
00337 {  
00338   Int_t xBins = 1;
00339   Double_t* xEdges = new Double_t[2];  
00340 
00341   xEdges[0] = -1.0e6;
00342   xEdges[1] = +1.0e6;
00343 
00344   Binning* newBinning = new Binning("",xBins,xEdges);
00345 
00346   delete [] xEdges;
00347 
00348   return newBinning; 
00349 }
00350     
00351 Binning* BinningScheme::BuildBeamCC()
00352 {
00353   //  assert(0);
00354   Int_t xBins = 100;
00355   Double_t* xEdges = new Double_t[101];
00356 
00357   // Binning Scheme
00358   // ==============
00359   // 0.5  GeV bins, up to 0.5 GeV
00360   // 0.25 GeV bins, up to 20  GeV
00361   // 1.0  GeV bins, up to 30  GeV
00362   // 2.0  GeV bins, up to 50  GeV
00363   // 150  GeV bin,  up to 200 GeV
00364 
00365   xEdges[0]  =  0.0;
00366   xEdges[1]  =  0.5;
00367 
00368   xEdges[2]  =  0.75;
00369   xEdges[3]  =  1.00;
00370   xEdges[4]  =  1.25;
00371   xEdges[5]  =  1.50;
00372   xEdges[6]  =  1.75;
00373   xEdges[7]  =  2.00;
00374   xEdges[8]  =  2.25;
00375   xEdges[9]  =  2.50;
00376   xEdges[10] =  2.75;
00377   xEdges[11] =  3.00;
00378   xEdges[12] =  3.25;
00379   xEdges[13] =  3.50;
00380   xEdges[14] =  3.75;
00381   xEdges[15] =  4.00;
00382   xEdges[16] =  4.25;
00383   xEdges[17] =  4.50;
00384   xEdges[18] =  4.75;
00385   xEdges[19] =  5.00;
00386   xEdges[20] =  5.25;
00387   xEdges[21] =  5.50;
00388   xEdges[22] =  5.75;
00389   xEdges[23] =  6.00;
00390   xEdges[24] =  6.25;
00391   xEdges[25] =  6.50;
00392   xEdges[26] =  6.75;
00393   xEdges[27] =  7.00;
00394   xEdges[28] =  7.25;
00395   xEdges[29] =  7.50;
00396   xEdges[30] =  7.75;
00397   xEdges[31] =  8.00;
00398   xEdges[32] =  8.25;
00399   xEdges[33] =  8.50;
00400   xEdges[34] =  8.75;
00401   xEdges[35] =  9.00;
00402   xEdges[36] =  9.25;
00403   xEdges[37] =  9.50;
00404   xEdges[38] =  9.75;
00405   xEdges[39] = 10.00;
00406   xEdges[40] = 10.25;
00407   xEdges[41] = 10.50;
00408   xEdges[42] = 10.75;
00409   xEdges[43] = 11.00;
00410   xEdges[44] = 11.25;
00411   xEdges[45] = 11.50;
00412   xEdges[46] = 11.75;
00413   xEdges[47] = 12.00;
00414   xEdges[48] = 12.25;
00415   xEdges[49] = 12.50;
00416   xEdges[50] = 12.75;
00417   xEdges[51] = 13.00;
00418   xEdges[52] = 13.25;
00419   xEdges[53] = 13.50;
00420   xEdges[54] = 13.75;
00421   xEdges[55] = 14.00;
00422   xEdges[56] = 14.25;
00423   xEdges[57] = 14.50;
00424   xEdges[58] = 14.75;
00425   xEdges[59] = 15.00;
00426   xEdges[60] = 15.25;
00427   xEdges[61] = 15.50;
00428   xEdges[62] = 15.75;
00429   xEdges[63] = 16.00;
00430   xEdges[64] = 16.25;
00431   xEdges[65] = 16.50;
00432   xEdges[66] = 16.75;
00433   xEdges[67] = 17.00;
00434   xEdges[68] = 17.25;
00435   xEdges[69] = 17.50;
00436   xEdges[70] = 17.75;
00437   xEdges[71] = 18.00;
00438   xEdges[72] = 18.25;
00439   xEdges[73] = 18.50;
00440   xEdges[74] = 18.75;
00441   xEdges[75] = 19.00;
00442   xEdges[76] = 19.25;
00443   xEdges[77] = 19.50;
00444   xEdges[78] = 19.75;
00445   xEdges[79] = 20.00;
00446 
00447   xEdges[80] = 21.0;
00448   xEdges[81] = 22.0;
00449   xEdges[82] = 23.0;
00450   xEdges[83] = 24.0;
00451   xEdges[84] = 25.0;
00452   xEdges[85] = 26.0;
00453   xEdges[86] = 27.0;
00454   xEdges[87] = 28.0;
00455   xEdges[88] = 29.0;
00456   xEdges[89] = 30.0;
00457 
00458   xEdges[90] = 32.0;
00459   xEdges[91] = 34.0;
00460   xEdges[92] = 36.0;
00461   xEdges[93] = 38.0;
00462   xEdges[94] = 40.0;
00463   xEdges[95] = 42.0;
00464   xEdges[96] = 44.0;
00465   xEdges[97] = 46.0;
00466   xEdges[98] = 48.0;
00467   xEdges[99] = 50.0;
00468 
00469   xEdges[100] = 200.0;
00470 
00471   Binning* newBinning = new Binning("Reconstructed E_{#nu} (GeV)",xBins,xEdges);
00472 
00473   delete [] xEdges;
00474 
00475   return newBinning;  
00476 }
00477 
00478 Binning* BinningScheme::BuildBeamCCMinosPlus()
00479 {
00480   Int_t xBins = 98;
00481   // Binning Scheme
00482   // ==============
00483   // 1.0  GeV bins, up to 1 GeV
00484   // 0.25 GeV bins, up to 20  GeV
00485   // 1.0  GeV bins, up to 30  GeV
00486   // 2.0  GeV bins, up to 50  GeV
00487   // 150  GeV bin,  up to 200 GeV
00488   std::vector<Double_t> edges;
00489   edges.push_back(0.0);
00490   Double_t E = 0;
00491   while (E<200) {
00492     if (E < 0.01) E+=1.0;
00493     else if (E < 19.99) E+=0.25;
00494     else if (E < 29.99) E+=1.0;
00495     else if (E < 49.99) E+=2.0;
00496     else if (E < 199.99) E+=150.0;
00497     edges.push_back(E);
00498   }
00499 
00500   assert(edges.size() == xBins+1);
00501   Binning* newBinning = new Binning("Reconstructed E_{#nu} (GeV)",xBins,&(edges.at(0))); //convert vector to c-style array by pointing at 1st element
00502   return newBinning;  
00503 }
00504 
00505 Binning* BinningScheme::BuildBeamNUE()
00506 {  
00507   Int_t xBins = 5;
00508   Double_t* xEdges = new Double_t[6];
00509 
00510   xEdges[0] = 1.0;
00511   xEdges[1] = 2.0;  
00512   xEdges[2] = 3.0;
00513   xEdges[3] = 4.0;
00514   xEdges[4] = 5.0;
00515   xEdges[5] = 8.0;
00516 
00517   Binning* newBinning = new Binning("Reconstructed E_{#nu} (GeV)",xBins,xEdges);
00518 
00519   delete [] xEdges;
00520 
00521   return newBinning;  
00522 }
00523     
00524 Binning* BinningScheme::BuildBeamRAF()
00525 {
00526   Int_t xBins = 123;
00527   Double_t* xEdges = new Double_t[124];
00528 
00529   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00530     xEdges[ibin] = (Double_t)ibin;
00531   }
00532 
00533   Binning* newBinning = new Binning("RAF Bin Number",xBins,xEdges);
00534 
00535   delete [] xEdges;
00536 
00537   return newBinning;  
00538 }
00539 
00540 Binning* BinningScheme::BuildBeamHiRes()
00541 {
00542   Int_t xBins = 120;
00543   Double_t* xEdges = new Double_t[121];
00544 
00545   // Binning Scheme
00546   // ==============
00547   // 0.5   GeV bins, up to   0.5 GeV
00548   // 0.25  GeV bins, up to   1.0 GeV
00549   // 0.125 GeV bins, up to   6.0 GeV
00550   // 0.25  GeV bins, up to  20   GeV
00551   // 1.0   GeV bins, up to  30   GeV
00552   // 2.0   GeV bins, up to  50   GeV
00553   // 150   GeV bin,  up to 200   GeV
00554 
00555   xEdges[0]  =  0.0;
00556   xEdges[1]  =  0.5;
00557 
00558   xEdges[2]   = 0.75;
00559   xEdges[3]   = 1.00;
00560 
00561   xEdges[4]   = 1.125;
00562   xEdges[5]   = 1.250;
00563   xEdges[6]   = 1.375;
00564   xEdges[7]   = 1.500;
00565   xEdges[8]   = 1.625;
00566   xEdges[9]   = 1.750;
00567   xEdges[10]  = 1.875;
00568   xEdges[11]  = 2.000;
00569   xEdges[12]  = 2.125;
00570   xEdges[13]  = 2.250;
00571   xEdges[14]  = 2.375;
00572   xEdges[15]  = 2.500;
00573   xEdges[16]  = 2.625;
00574   xEdges[17]  = 2.750;
00575   xEdges[18]  = 2.875;
00576   xEdges[19]  = 3.000;
00577   xEdges[20]  = 3.125;
00578   xEdges[21]  = 3.250;
00579   xEdges[22]  = 3.375;
00580   xEdges[23]  = 3.500;
00581   xEdges[24]  = 3.625;
00582   xEdges[25]  = 3.750;
00583   xEdges[26]  = 3.875;
00584   xEdges[27]  = 4.000;
00585   xEdges[28]  = 4.125;
00586   xEdges[29]  = 4.250;
00587   xEdges[30]  = 4.375;
00588   xEdges[31]  = 4.500;
00589   xEdges[32]  = 4.625;
00590   xEdges[33]  = 4.750;
00591   xEdges[34]  = 4.875;
00592   xEdges[35]  = 5.000;
00593   xEdges[36]  = 5.125;
00594   xEdges[37]  = 5.250;
00595   xEdges[38]  = 5.375;
00596   xEdges[39]  = 5.500;
00597   xEdges[40]  = 5.625;
00598   xEdges[41]  = 5.750;
00599   xEdges[42]  = 5.875;
00600   xEdges[43]  = 6.000;
00601 
00602   xEdges[44] =  6.25;
00603   xEdges[45] =  6.50;
00604   xEdges[46] =  6.75;
00605   xEdges[47] =  7.00;
00606   xEdges[48] =  7.25;
00607   xEdges[49] =  7.50;
00608   xEdges[50] =  7.75;
00609   xEdges[51] =  8.00;
00610   xEdges[52] =  8.25;
00611   xEdges[53] =  8.50;
00612   xEdges[54] =  8.75;
00613   xEdges[55] =  9.00;
00614   xEdges[56] =  9.25;
00615   xEdges[57] =  9.50;
00616   xEdges[58] =  9.75;
00617   xEdges[59] = 10.00;
00618   xEdges[60] = 10.25;
00619   xEdges[61] = 10.50;
00620   xEdges[62] = 10.75;
00621   xEdges[63] = 11.00;
00622   xEdges[64] = 11.25;
00623   xEdges[65] = 11.50;
00624   xEdges[66] = 11.75;
00625   xEdges[67] = 12.00;
00626   xEdges[68] = 12.25;
00627   xEdges[69] = 12.50;
00628   xEdges[70] = 12.75;
00629   xEdges[71] = 13.00;
00630   xEdges[72] = 13.25;
00631   xEdges[73] = 13.50;
00632   xEdges[74] = 13.75;
00633   xEdges[75] = 14.00;
00634   xEdges[76] = 14.25;
00635   xEdges[77] = 14.50;
00636   xEdges[78] = 14.75;
00637   xEdges[79] = 15.00;
00638   xEdges[80] = 15.25;
00639   xEdges[81] = 15.50;
00640   xEdges[82] = 15.75;
00641   xEdges[83] = 16.00;
00642   xEdges[84] = 16.25;
00643   xEdges[85] = 16.50;
00644   xEdges[86] = 16.75;
00645   xEdges[87] = 17.00;
00646   xEdges[88] = 17.25;
00647   xEdges[89] = 17.50;
00648   xEdges[90] = 17.75;
00649   xEdges[91] = 18.00;
00650   xEdges[92] = 18.25;
00651   xEdges[93] = 18.50;
00652   xEdges[94] = 18.75;
00653   xEdges[95] = 19.00;
00654   xEdges[96] = 19.25;
00655   xEdges[97] = 19.50;
00656   xEdges[98] = 19.75;
00657   xEdges[99] = 20.00;
00658 
00659   xEdges[100] = 21.0;
00660   xEdges[101] = 22.0;
00661   xEdges[102] = 23.0;
00662   xEdges[103] = 24.0;
00663   xEdges[104] = 25.0;
00664   xEdges[105] = 26.0;
00665   xEdges[106] = 27.0;
00666   xEdges[107] = 28.0;
00667   xEdges[108] = 29.0;
00668   xEdges[109] = 30.0;
00669 
00670   xEdges[110] = 32.0;
00671   xEdges[111] = 34.0;
00672   xEdges[112] = 36.0;
00673   xEdges[113] = 38.0;
00674   xEdges[114] = 40.0;
00675   xEdges[115] = 42.0;
00676   xEdges[116] = 44.0;
00677   xEdges[117] = 46.0;
00678   xEdges[118] = 48.0;
00679   xEdges[119] = 50.0;
00680 
00681   xEdges[120] = 200.0;
00682 
00683   Binning* newBinning = new Binning("Reconstructed E_{#nu} (GeV)",xBins,xEdges);
00684 
00685   delete [] xEdges;
00686 
00687   return newBinning;  
00688 }
00689 
00690 Binning* BinningScheme::BuildAtmosSingleBin()
00691 {
00692   Int_t xBins = 1;
00693   Double_t* xEdges = new Double_t[2];  
00694 
00695   xEdges[0] = pow(10.0,-0.3); // matching
00696   xEdges[1] = pow(10.0,+1.3); // kAtmosEnergy
00697 
00698   Binning* newBinning = new Binning("",xBins,xEdges);
00699 
00700   delete [] xEdges;
00701 
00702   return newBinning; 
00703 }
00704 
00705 Binning* BinningScheme::BuildAtmosLogLE()
00706 {
00707   Int_t xBins = 25;
00708   Double_t* xEdges = new Double_t[26];
00709 
00710   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00711     xEdges[ibin] = pow(10.0,-0.5+5.0*ibin/double(xBins));
00712     // std::cout << " --- debug: [" << ibin << "] " << -0.5+5.0*ibin/double(xBins) << std::endl;
00713   }
00714 
00715   Binning* newBinning = new Binning("Reconstructed L_{#nu}/E_{#nu} (km/GeV)",xBins,xEdges);
00716 
00717   delete [] xEdges;
00718 
00719   return newBinning;  
00720 }
00721    
00722 Binning* BinningScheme::BuildAtmosEnergy()
00723 {
00724   Int_t xBins = 16;
00725   Double_t* xEdges = new Double_t[17];
00726 
00727   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00728     xEdges[ibin] = pow(10.0,-0.3+1.6*ibin/double(xBins));
00729   }
00730 
00731   Binning* newBinning = new Binning("Reconstructed E_{#nu} (GeV)",xBins,xEdges);
00732 
00733   delete [] xEdges;
00734 
00735   return newBinning;  
00736 } 
00737 
00738 Binning* BinningScheme::BuildAtmosAngle()
00739 {
00740   Int_t xBins = 30;
00741   Double_t* xEdges = new Double_t[31];
00742 
00743   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00744     xEdges[ibin] = -1.0 + 2.0*ibin/double(xBins);
00745   }
00746 
00747   Binning* newBinning = new Binning("cos(#theta_{z})",xBins,xEdges);
00748 
00749   delete [] xEdges;
00750 
00751   return newBinning;  
00752 }
00753 
00754 Binning* BinningScheme::BuildAtmosAngleLoRes()
00755 {
00756   Int_t xBins = 20;
00757   Double_t* xEdges = new Double_t[21];
00758 
00759   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00760     xEdges[ibin] = -1.0 + 2.0*ibin/double(xBins);
00761   }
00762 
00763   Binning* newBinning = new Binning("cos(#theta_{z})",xBins,xEdges);
00764 
00765   delete [] xEdges;
00766 
00767   return newBinning;  
00768 }
00769 
00770 Binning* BinningScheme::BuildAtmosAngleHiRes()
00771 {
00772   Int_t xBins = 40;
00773   Double_t* xEdges = new Double_t[41];
00774 
00775   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00776     xEdges[ibin] = -1.0 + 2.0*ibin/double(xBins);
00777   }
00778 
00779   Binning* newBinning = new Binning("cos(#theta_{z})",xBins,xEdges);
00780 
00781   delete [] xEdges;
00782 
00783   return newBinning;  
00784 }
00785     
00786 Binning* BinningScheme::BuildAtmos2D()
00787 {
00788   Int_t xBins = 30;
00789   Double_t* xEdges = new Double_t[31];
00790 
00791   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00792     xEdges[ibin] = -1.0 + 2.0*ibin/double(xBins);
00793   }
00794 
00795   Int_t yBins = 36;
00796   Double_t* yEdges = new Double_t[37];
00797 
00798   for( Int_t ibin=0; ibin<yBins+1; ibin++ ){
00799     yEdges[ibin] = pow(10.0,-0.6+2.4*ibin/double(yBins));
00800   }
00801 
00802   Binning* newBinning = new Binning("cos(#theta_{z})",xBins,xEdges,
00803                                     "Reconstructed E_{#nu} (GeV)",yBins,yEdges);
00804   delete [] xEdges;
00805   delete [] yEdges;
00806 
00807   return newBinning;  
00808 }
00809 
00810 Binning* BinningScheme::BuildAtmos2DLoRes()
00811 {
00812   Int_t xBins = 20;
00813   Double_t* xEdges = new Double_t[21];
00814 
00815   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00816     xEdges[ibin] = -1.0 + 2.0*ibin/double(xBins);
00817   }
00818 
00819   Int_t yBins = 24;
00820   Double_t* yEdges = new Double_t[25];
00821 
00822   for( Int_t ibin=0; ibin<yBins+1; ibin++ ){
00823     yEdges[ibin] = pow(10.0,-0.6+2.4*ibin/double(yBins));
00824   }
00825 
00826   Binning* newBinning = new Binning("cos(#theta_{z})",xBins,xEdges,
00827                                     "Reconstructed E_{#nu} (GeV)",yBins,yEdges);
00828   delete [] xEdges;
00829   delete [] yEdges;
00830 
00831   return newBinning;  
00832 }
00833 
00834 Binning* BinningScheme::BuildAtmos2DHiRes()
00835 {
00836   Int_t xBins = 40;
00837   Double_t* xEdges = new Double_t[41];
00838 
00839   for( Int_t ibin=0; ibin<xBins+1; ibin++ ){
00840     xEdges[ibin] = -1.0 + 2.0*ibin/double(xBins);
00841   }
00842 
00843   /*
00844   Int_t yBins = 48;
00845   Double_t* yEdges = new Double_t[49];
00846 
00847   for( Int_t ibin=0; ibin<yBins+1; ibin++ ){
00848     yEdges[ibin] = pow(10.0,-0.6+2.4*ibin/double(yBins));
00849   }
00850   */
00851 
00852   Int_t yBins = 60;
00853   Double_t* yEdges = new Double_t[61];
00854 
00855   for( Int_t ibin=0; ibin<yBins+1; ibin++ ){
00856     yEdges[ibin] = pow(10.0,-1.0+3.0*ibin/double(yBins));
00857   }
00858 
00859   Binning* newBinning = new Binning("cos(#theta_{z})",xBins,xEdges,
00860                                     "Reconstructed E_{#nu} (GeV)",yBins,yEdges);
00861   delete [] xEdges;
00862   delete [] yEdges;
00863 
00864   return newBinning;  
00865 }
00866 
00867 void BinningScheme::ReadBinningScheme( const char* filename )
00868 {
00869   std::cout << " *** BinningScheme::ReadBinningScheme(...) *** " << std::endl
00870             << "  filename = " << filename << std::endl;
00871 
00872   if( fBinningSchemeIsLoaded ) {
00873     std::cout << "  warning: binning scheme is already loaded " << std::endl;
00874     return;
00875   }
00876 
00877   if( gSystem->AccessPathName(filename) ) { // returns kTrue if non-existent
00878     std::cout << "  warning: failed to find " << filename << std::endl;
00879     return;
00880   }
00881 
00882   BinningScheme_t myScheme = kUnknownScheme;
00883 
00884   TDirectory* tmpd = gDirectory;
00885   TFile* file = new TFile(filename,"READ");
00886 
00887   TH1D* hBinning = (TH1D*)(file->Get("hBinning"));
00888 
00889   if( hBinning ){
00890     if( hBinning->GetBinContent(1)>-0.5 
00891      && hBinning->GetBinContent(1)<+0.5 )   myScheme = kScheme0;
00892     if( hBinning->GetBinContent(1)>+0.5 
00893      && hBinning->GetBinContent(1)<+1.5 )   myScheme = kScheme1;
00894     if( hBinning->GetBinContent(1)>+1.5 
00895      && hBinning->GetBinContent(1)<+2.5 )   myScheme = kScheme2;
00896     if( hBinning->GetBinContent(1)>+99.5 
00897      && hBinning->GetBinContent(1)<+100.5 ) myScheme = kCounterForMINOS;
00898     if( hBinning->GetBinContent(1)>+100.5 
00899      && hBinning->GetBinContent(1)<+101.5 ) myScheme = kCounterForLBNE;
00900   }
00901 
00902   file->Close();
00903   tmpd->cd();
00904 
00905   assert( myScheme!= kUnknownScheme );
00906 
00907   std::cout << "  binning scheme: " << AsString(myScheme) << std::endl;
00908 
00909   SetBinningScheme( myScheme );
00910 
00911   return;
00912 }
00913     
00914 void BinningScheme::WriteBinningScheme( const char* filename )
00915 {
00916   std::cout << " *** BinningScheme::WriteBinningScheme(...) *** " << std::endl
00917             << "  filename = " << filename << std::endl;
00918 
00919   TH1D* hBinning = new TH1D("hBinning","",1,0,1);
00920 
00921   if( fBinningScheme==kScheme0         ) hBinning->SetBinContent(1,0);
00922   if( fBinningScheme==kScheme1         ) hBinning->SetBinContent(1,1);
00923   if( fBinningScheme==kScheme1A        ) hBinning->SetBinContent(1,1);
00924   if( fBinningScheme==kScheme2         ) hBinning->SetBinContent(1,2);
00925   if( fBinningScheme==kScheme2A        ) hBinning->SetBinContent(1,2);
00926   if( fBinningScheme==kCounterForMINOS ) hBinning->SetBinContent(1,100);
00927   if( fBinningScheme==kCounterForLBNE  ) hBinning->SetBinContent(1,101);
00928  
00929   std::cout << "  binning scheme = " << hBinning->GetBinContent(1) << std::endl;
00930 
00931   TDirectory* tmpd = gDirectory;
00932   TFile* file = new TFile(filename,"RECREATE");
00933 
00934   hBinning->SetDirectory(file);
00935   hBinning->Write();
00936 
00937   file->Close();
00938   tmpd->cd();
00939 
00940   std::cout << "  Done! " << std::endl;
00941 
00942   return;
00943 }