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
00083
00084
00085
00086 if( run==kAtmos ){
00087 if( resolution==kLoRes ){
00088 return kSingleBin;
00089 }
00090 else{
00091 return kAtmosLogLE;
00092 }
00093 }
00094
00095
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
00117
00118
00119
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
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
00165
00166
00167
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
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;
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
00354 Int_t xBins = 100;
00355 Double_t* xEdges = new Double_t[101];
00356
00357
00358
00359
00360
00361
00362
00363
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
00482
00483
00484
00485
00486
00487
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)));
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
00546
00547
00548
00549
00550
00551
00552
00553
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);
00696 xEdges[1] = pow(10.0,+1.3);
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
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
00845
00846
00847
00848
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) ) {
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 }