AtNuResolution Class Reference

#include <AtNuResolution.h>

List of all members.

Public Member Functions

 AtNuResolution ()
 ~AtNuResolution ()
Double_t GetFluxWeight (Int_t fIDnu, Double_t fEtrue, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior=AtNuResolutionPrior::kFlat)
Double_t GetNeugenWeight (Int_t fIDnu, Double_t fEtrue, Double_t fYtrue, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior=AtNuResolutionPrior::kFlat)
Double_t GetResolution (AtNuResolutionEventType::AtNuResolutionEventType_t eventtype, Bool_t gooddirection, Bool_t goodenergy, Bool_t rangecurve, Int_t fCharge=0, Double_t fCosX=0.0, Double_t fCosY=0.0, Double_t fCosZ=0.0, Double_t fEmu=0.0, Double_t fEmuErr=0.0, Double_t fEshw=0.0, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior=AtNuResolutionPrior::kFlat)
TH1D * GetResolutionPdf (AtNuResolutionEventType::AtNuResolutionEventType_t eventtype, Bool_t gooddirection, Bool_t goodenergy, Bool_t rangecurve, Int_t fCharge=0, Double_t fCosX=0.0, Double_t fCosY=0.0, Double_t fCosZ=0.0, Double_t fEmu=0.0, Double_t fEmuErr=0.0, Double_t fEshw=0.0, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior=AtNuResolutionPrior::kFlat, Bool_t renormalise=0)
AtNuResolutionType::AtNuResolutionType_t GetResolutionType (AtNuResolutionEventType::AtNuResolutionEventType_t eventtype, Bool_t gooddirection, Bool_t goodenergy, Bool_t rangecurve, Int_t fCharge=0, Double_t fCosX=0.0, Double_t fCosY=0.0, Double_t fCosZ=0.0, Double_t fEmu=0.0, Double_t fEmuErr=0.0, Double_t fEshw=0.0, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior=AtNuResolutionPrior::kFlat)
AtNuResolutionType::AtNuResolutionType_t GetResolutionType (AtNuResolutionEventType::AtNuResolutionEventType_t eventtype, Double_t resolution)
void SetEnergyBinning (Int_t bins, Double_t min, Double_t max)
void LoadHistograms (const char *filename, Bool_t flag=1)

Private Member Functions

Double_t GetResolutionCV (Bool_t gooddirection, Bool_t goodenergy, Bool_t rangecurve, Int_t fCharge=0, Double_t fCosX=0.0, Double_t fCosY=0.0, Double_t fCosZ=0.0, Double_t fEmu=0.0, Double_t fEmuErr=0.0, Double_t fEshw=0.0, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior=AtNuResolutionPrior::kFlat)
Double_t GetResolutionUPMU (Bool_t gooddirection, Bool_t goodenergy, Bool_t rangecurve, Int_t fCharge=0, Double_t fCosX=0.0, Double_t fCosY=0.0, Double_t fCosZ=0.0, Double_t fEmu=0.0, Double_t fEmuErr=0.0, Double_t fEshw=0.0, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior=AtNuResolutionPrior::kFlat)
void AddToBayesPdf (AtNuResBayesPdf *&bayespdf, AtNuResNeugenPdf *neugenpdf, Double_t fEtrueWeight, Bool_t rangecurve, Double_t fEtrue, Double_t fCosX, Double_t fCosY, Double_t fCosZ, Double_t fEmu, Double_t fEmuErr, Double_t fEshw, Bool_t isQE, AtNuResolutionPrior::AtNuResolutionPrior_t bayesprior)
void LoadHistogramsFromCambridge (const char *filename)

Private Attributes

AtNuResBayesPdffBayesPdf
AtNuResEnergyPdffEnergyPdf
AtNuResNeugenPdffNeugenPdf
AtNuResNeugenPdf ** fNuQE
AtNuResNeugenPdf ** fNuRES
AtNuResNeugenPdf ** fNuDIS
AtNuResNeugenPdf ** fNuBarQE
AtNuResNeugenPdf ** fNuBarRES
AtNuResNeugenPdf ** fNuBarDIS
Double_t * fEnergyWeight
Double_t * fEnergyMedian
Bool_t fUseBayes
Int_t fHistNumber

Detailed Description

Definition at line 109 of file AtNuResolution.h.


Constructor & Destructor Documentation

AtNuResolution::AtNuResolution (  ) 

Definition at line 12 of file AtNuResolution.cxx.

References n.

00013 {
00014   std::cout << " *** AtNuResolution::AtNuResolution() *** " << std::endl;
00015 
00016   // neutrino weights
00017   // ================
00018   fNuQE = new AtNuResNeugenPdf*[1000];
00019   for( Int_t n=0; n<1000; n++ ){
00020     fNuQE[n] = 0;
00021   }
00022 
00023   fNuRES = new AtNuResNeugenPdf*[1000];
00024   for( Int_t n=0; n<1000; n++ ){
00025     fNuRES[n] = 0;
00026   }
00027 
00028   fNuDIS = new AtNuResNeugenPdf*[1000];
00029   for( Int_t n=0; n<1000; n++ ){
00030     fNuDIS[n] = 0;
00031   }
00032 
00033   // anti-neutrino weights
00034   // =====================
00035   fNuBarQE = new AtNuResNeugenPdf*[1000];
00036   for( Int_t n=0; n<1000; n++ ){
00037     fNuBarQE[n] = 0;
00038   }
00039 
00040   fNuBarRES = new AtNuResNeugenPdf*[1000];
00041   for( Int_t n=0; n<1000; n++ ){
00042     fNuBarRES[n] = 0;
00043   }
00044 
00045   fNuBarDIS = new AtNuResNeugenPdf*[1000];
00046   for( Int_t n=0; n<1000; n++ ){
00047     fNuBarDIS[n] = 0;
00048   }
00049 
00050   // energy weighting
00051   // ================
00052   fEnergyWeight = new Double_t[1000];
00053   for( Int_t n=0; n<1000; n++ ){
00054     fEnergyWeight[n] = 0.0;
00055   }
00056 
00057   fEnergyMedian = new Double_t[1000];
00058   for( Int_t n=0; n<1000; n++ ){
00059     fEnergyMedian[n] = 0.0;
00060   }
00061 
00062   // bayes pdf
00063   // =========
00064   fHistNumber = 0;
00065   fBayesPdf = 0;
00066 
00067   // energy pdfs
00068   // ===========
00069   fEnergyPdf = new AtNuResEnergyPdf();
00070 
00071   // neugen pdfs (loaded)
00072   // ====================
00073   fUseBayes = 0;
00074 }

AtNuResolution::~AtNuResolution (  ) 

Definition at line 76 of file AtNuResolution.cxx.

References fBayesPdf, fEnergyMedian, fEnergyPdf, fEnergyWeight, fNuBarDIS, fNuBarQE, fNuBarRES, fNuDIS, fNuQE, fNuRES, and n.

00077 {
00078   std::cout << " *** AtNuResolution::~AtNuResolution() *** " << std::endl;
00079 
00080   // neutrino weights
00081   // ================
00082   for( Int_t n=0; n<1000; n++ ){
00083     if( fNuQE[n] ) delete fNuQE[n];
00084   }
00085   delete [] fNuQE;
00086 
00087   for( Int_t n=0; n<1000; n++ ){
00088     if( fNuRES[n] ) delete fNuRES[n];
00089   }
00090   delete [] fNuRES;
00091 
00092   for( Int_t n=0; n<1000; n++ ){
00093     if( fNuDIS[n] ) delete fNuDIS[n];
00094   }
00095   delete [] fNuDIS;
00096 
00097   // anti-neutrino weights
00098   // =====================
00099   for( Int_t n=0; n<1000; n++ ){
00100     if( fNuBarQE[n] ) delete fNuBarQE[n];
00101   }
00102   delete [] fNuBarQE;
00103 
00104   for( Int_t n=0; n<1000; n++ ){
00105     if( fNuBarRES[n] ) delete fNuBarRES[n];
00106   }
00107   delete [] fNuBarRES;
00108 
00109   for( Int_t n=0; n<1000; n++ ){
00110     if( fNuBarDIS[n] ) delete fNuBarDIS[n];
00111   }
00112   delete [] fNuBarDIS;
00113 
00114   // energy weighting
00115   // ================
00116   delete [] fEnergyWeight;
00117   delete [] fEnergyMedian;
00118 
00119   // energy pdf
00120   // ==========
00121   if( fEnergyPdf ) delete fEnergyPdf;
00122 
00123   // bayes pdf
00124   // =========
00125   if( fBayesPdf ) delete fBayesPdf;
00126 }  


Member Function Documentation

void AtNuResolution::AddToBayesPdf ( AtNuResBayesPdf *&  bayespdf,
AtNuResNeugenPdf neugenpdf,
Double_t  fEtrueWeight,
Bool_t  rangecurve,
Double_t  fEtrue,
Double_t  fCosX,
Double_t  fCosY,
Double_t  fCosZ,
Double_t  fEmu,
Double_t  fEmuErr,
Double_t  fEshw,
Bool_t  isQE,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior 
) [private]

Definition at line 448 of file AtNuResolution.cxx.

References fEnergyPdf, AtNuResBayesPdf::Fill(), GetFluxWeight(), AtNuResNeugenPdf::GetIDnu(), AtNuResEnergyPdf::GetMuonCurveWeightCorrected(), AtNuResEnergyPdf::GetMuonRangeWeightCorrected(), GetNeugenWeight(), AtNuResEnergyPdf::GetShowerWeightCorrected(), AtNuResNeugenPdf::GetW2(), AtNuResNeugenPdf::GetW2bins(), AtNuResNeugenPdf::GetWeight(), AtNuResNeugenPdf::GetY(), AtNuResNeugenPdf::GetYbins(), and Mphysical::pi.

Referenced by GetResolutionCV().

00449 {
00450   // sanity check
00451   if( neugenpdf==0 ){
00452     std::cout << " *** Warning: Failed to load Neugen pdf *** " << std::endl;
00453     return;
00454   }
00455 
00456   Double_t Y,W2;
00457   Double_t Emu,Eshw;
00458   Double_t Pmu,Pshw;
00459   Double_t CosTH,SinTH;
00460   Double_t CosPH,SinPH;
00461   Double_t fCosTH,fSinTH;
00462   Double_t fCosPH,fSinPH;
00463   Double_t muCosX,muCosY,muCosZ;
00464   Double_t nuCosX,nuCosY,nuCosZ;
00465   Double_t nuE,nuL,logLE;
00466 
00467   Double_t weight;
00468   Double_t priorweight;
00469   Double_t muonweight;
00470   Double_t showerweight;
00471   Double_t neugenweight;
00472   Double_t solidangleweight;
00473 
00474   Int_t fPhiBins = 60;
00475   Double_t muMass = 0.106;     // muon mass (GeV)
00476   Double_t pMass = 0.938;      // proton mass (GeV)
00477   Double_t pi = 3.1415926535;  // pi
00478   
00479   Double_t fR = 6371.0;        // radius of earth
00480   Double_t fH = 15.0;          // production height
00481   Double_t fD = 0.65;          // detector depth
00482 
00483   Double_t goodneugenweight = 0.0;
00484   Double_t totalneugenweight = 0.0;
00485 
00486   // neutrino type
00487   Int_t fIDnu = neugenpdf->GetIDnu();
00488 
00489   // observed spherical coordinates
00490   fCosTH = fCosY;
00491   fSinTH = sqrt(fCosX*fCosX+fCosZ*fCosZ);
00492   if( fSinTH>0.0 ){
00493     fCosPH = fCosX/fSinTH; 
00494     fSinPH = fCosZ/fSinTH;
00495   }
00496   else{
00497     fCosPH = 1.0; fSinPH = 0.0;
00498   }
00499 
00500   // loop over possible Y and W2 values
00501   for( Int_t nY=0; nY<neugenpdf->GetYbins(); nY++ ){
00502     for( Int_t nW2=0; nW2<neugenpdf->GetW2bins(); nW2++ ){
00503 
00504       // weight for neutrino interaction kinematics
00505       Y = neugenpdf->GetY(nY);
00506       W2 = neugenpdf->GetW2(nW2);
00507       neugenweight = neugenpdf->GetWeight(nY,nW2);
00508       totalneugenweight += neugenweight;
00509 
00510       // weight for flux and efficiency
00511       priorweight = fEtrueWeight
00512                   * this->GetFluxWeight(fIDnu,fEtrue,bayesprior)
00513                   * this->GetNeugenWeight(fIDnu,fEtrue,Y,bayesprior);
00514  
00515       // extract muon and shower kinematics
00516       Eshw = Y*fEtrue+pMass;
00517       Emu = (1.0-Y)*fEtrue;
00518 
00519       // require physical muon and shower energies
00520       if( Eshw*Eshw>W2 && Emu>muMass ){
00521 
00522         Pshw = sqrt(Eshw*Eshw-W2);
00523         Pmu  = sqrt(Emu*Emu-muMass*muMass);
00524 
00525         // require physical muon and shower momenta
00526         if( fEtrue>fabs(Pmu-Pshw) 
00527          && fEtrue<fabs(Pmu+Pshw) ){
00528           goodneugenweight += neugenweight;
00529 
00530           // weight for muon energy resolution 
00531           if( rangecurve ){
00532             muonweight = fEnergyPdf->GetMuonRangeWeightCorrected(fEmu,Emu); // (ereco,etrue)
00533           }
00534           else{
00535             muonweight = fEnergyPdf->GetMuonCurveWeightCorrected(fEmu,Emu,fEmuErr); // (ereco,etrue,error)
00536           }
00537  
00538           // weight for shower energy resolution
00539           showerweight = fEnergyPdf->GetShowerWeightCorrected(fEshw,fEtrue-Emu,fIDnu,isQE);  // (ereco,etrue,inu,isqe)
00540         
00541           // integrate over spherical coordinates
00542           CosTH = (fEtrue*fEtrue+Pmu*Pmu-Pshw*Pshw)/(2.0*Pmu*fEtrue);
00543           SinTH = sqrt(1.0-CosTH*CosTH);         
00544 
00545           // weight for solid angle
00546           solidangleweight = SinTH;
00547 
00548           // total weight (interaction and resolution)
00549           weight = priorweight*neugenweight*muonweight*showerweight*solidangleweight;
00550 
00551           // loop over polar angles
00552           for( Int_t nPhi=0; nPhi<fPhiBins; nPhi++ ){
00553             CosPH = cos(2.0*pi*nPhi/fPhiBins);
00554             SinPH = sin(2.0*pi*nPhi/fPhiBins);
00555 
00556             muCosX = SinTH*CosPH;
00557             muCosY = CosTH;
00558             muCosZ = SinTH*SinPH;
00559 
00560             nuCosX =   muCosX*fCosTH*fCosPH + muCosY*fSinTH*fCosPH - muCosZ*fSinPH;
00561             nuCosY = - muCosX*fSinTH        + muCosY*fCosTH;
00562             nuCosZ =   muCosX*fCosTH*fSinPH + muCosY*fSinTH*fSinPH + muCosZ*fCosPH;
00563 
00564             nuE = fEtrue;
00565             nuL = + (fR-fD)*nuCosY
00566                   + sqrt( (fR-fD)*(fR-fD)*nuCosY*nuCosY+(fH+fD)*(2.0*fR+fH-fD) );
00567 
00568             logLE = log10(nuL/nuE);
00569 
00570             //
00571             // std::cout << " neugenweight=" << neugenweight << " muonweight=" << muonweight << " showerweight=" << showerweight << std::endl;
00572             // std::cout << " nuCosY=" << nuCosY << " nuL=" << nuL << " nuE=" << nuE << " LogLE=" << logLE << " weight=" << weight << std::endl;
00573             //
00574 
00575             bayespdf->Fill(logLE,weight);  
00576 
00577           } // for( Int_t nPhi=0; nPhi<fPhiBins; nPhi++ ){...
00578         } // if( fEtrue>(Pmu+Pshw) ){...
00579       } // if( Eshw*Eshw>W2 && Emu>muMass ){...
00580 
00581     } // for( Int_t nW2=0; nW2<fW2bins; nW2++ ){...
00582   } // for( Int_t nY=0; nY<fYbins; nY++ ){...
00583   
00584   return;
00585 }

Double_t AtNuResolution::GetFluxWeight ( Int_t  fIDnu,
Double_t  fEtrue,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior = AtNuResolutionPrior::kFlat 
)

Definition at line 399 of file AtNuResolution.cxx.

References AtNuResolutionPrior::kAtmos.

Referenced by AddToBayesPdf().

00400 {
00401   // this method should return a weight that reflects
00402   // the relative rate of selected events at Etrue
00403     
00404   // sanity check
00405   if( fEtrue<=0 ){
00406     return 0.0;
00407   }
00408 
00409   // start with unity
00410   Double_t weight = 1.0;
00411   Double_t etrue = 0.0;
00412   
00413   // realistic spectrum
00414   // (a) use fermi function for rising edge of distribution
00415   // (b) use power law function for falling edge of distribution
00416   if( bayesprior==AtNuResolutionPrior::kAtmos ){
00417 
00418     if( fIDnu==+14 ){  // neutrinos (well-measured charge sign)  
00419       if( fEtrue>0.0 ){
00420         weight *= 1.0/(1.0+exp(-(fEtrue-0.90)/0.090));
00421       }
00422       if( fEtrue>1.5 ){
00423         etrue = fEtrue-(fEtrue-1.5)*exp(-(fEtrue-1.5)/0.25);
00424         weight *= 1.387/(1.0+0.389*(etrue-0.5)*(etrue-0.5)); 
00425       }
00426     }
00427 
00428     else if( fIDnu==-14 ){  // anti-neutrinos (well-measured charge sign)
00429       if( fEtrue>0.0 ){
00430         weight *= 1.0/(1.0+exp(-(fEtrue-0.90)/0.090));
00431       }
00432       if( fEtrue>1.5 ){
00433         etrue = fEtrue-0.5*exp(-(fEtrue-1.5)/0.75);
00434         weight *= 1.0/(1.0+0.30*(etrue-1.0)*(etrue-1.0)); 
00435       }
00436     }
00437 
00438     else {  // combined nu/nubar (no cuts on direction/charge)
00439       if( fEtrue>0.0 ) weight *= 1.0/(1.0+exp(-(fEtrue-0.55)/0.055));
00440       if( fEtrue>1.0 ) weight *= 1.0/(1.0+pow(fEtrue-1.0,1.7));
00441     }
00442   }
00443 
00444   // return weight
00445   return weight;
00446 }

Double_t AtNuResolution::GetNeugenWeight ( Int_t  fIDnu,
Double_t  fEtrue,
Double_t  fYtrue,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior = AtNuResolutionPrior::kFlat 
)

Definition at line 385 of file AtNuResolution.cxx.

Referenced by AddToBayesPdf().

00386 {
00387   // this method should return a weight that reflects
00388   // the relative rate of selected events at Ytrue
00389 
00390   // sanity check
00391   if( fEtrue<=0 ){
00392     return 0.0;
00393   }
00394 
00395   // for now, return unity
00396   return 1.0;
00397 }

Double_t AtNuResolution::GetResolution ( AtNuResolutionEventType::AtNuResolutionEventType_t  eventtype,
Bool_t  gooddirection,
Bool_t  goodenergy,
Bool_t  rangecurve,
Int_t  fCharge = 0,
Double_t  fCosX = 0.0,
Double_t  fCosY = 0.0,
Double_t  fCosZ = 0.0,
Double_t  fEmu = 0.0,
Double_t  fEmuErr = 0.0,
Double_t  fEshw = 0.0,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior = AtNuResolutionPrior::kFlat 
)

Definition at line 128 of file AtNuResolution.cxx.

References GetResolutionCV(), GetResolutionUPMU(), AtNuResolutionEventType::kCV, AtNuResolutionEventType::kNUE, and AtNuResolutionEventType::kUPMU.

Referenced by GetResolutionPdf(), GetResolutionType(), and AtNuResNtuple::Run().

00129 {
00130   Double_t resolution = 999.9;
00131 
00132   switch( eventtype ){
00133     case AtNuResolutionEventType::kNUE: 
00134       break;
00135     case AtNuResolutionEventType::kCV:
00136       resolution = GetResolutionCV(gooddirection,goodenergy,rangecurve,
00137                                    fCharge,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,bayesprior);
00138       break;
00139     case AtNuResolutionEventType::kUPMU:
00140       resolution = GetResolutionUPMU(gooddirection,goodenergy,rangecurve,
00141                                      fCharge,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,bayesprior);
00142       break;
00143     default:
00144       break;
00145   }
00146 
00147   //
00148   // std::cout << " *** LogLE Resolution = " << resolution << " *** " << std::endl;
00149   //
00150 
00151   return resolution;
00152 }

Double_t AtNuResolution::GetResolutionCV ( Bool_t  gooddirection,
Bool_t  goodenergy,
Bool_t  rangecurve,
Int_t  fCharge = 0,
Double_t  fCosX = 0.0,
Double_t  fCosY = 0.0,
Double_t  fCosZ = 0.0,
Double_t  fEmu = 0.0,
Double_t  fEmuErr = 0.0,
Double_t  fEshw = 0.0,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior = AtNuResolutionPrior::kFlat 
) [private]

Definition at line 265 of file AtNuResolution.cxx.

References AddToBayesPdf(), fBayesPdf, fEnergyMedian, fEnergyPdf, fEnergyWeight, fNuBarDIS, fNuBarQE, fNuBarRES, fNuDIS, fNuQE, fNuRES, fUseBayes, AtNuResBayesPdf::GetRMS(), AtNuResBayesPdf::GetTotalWeight(), and AtNuResBayesPdf::Reset().

Referenced by GetResolution().

00266 { 
00267   // Bayesian Resolution for Contained Vertex Events
00268   // ===============================================
00269   //
00270   // std::cout << " *** AtNuResolution::GetResolutionCV(...) *** " << std::endl;
00271   // std::cout << "   gooddirection=" << gooddirection << ", goodenergy=" << goodenergy << ", rangecurve=" << rangecurve << std::endl << "   charge=" << fCharge << ", cosX=" << fCosX << ", cosY=" << fCosY << ", cosZ=" << fCosZ << ", emu=" << fEmu << ", emuerr=" << fEmuErr << ", eshw=" << fEshw << std::endl;
00272   //
00273 
00274   // Reset Bayes Pdf
00275   // ===============
00276   if( fBayesPdf ) fBayesPdf->Reset();
00277 
00278   // bad direction
00279   // =============
00280   if( !gooddirection ){
00281     return 999.9; // return low resolution
00282   }
00283 
00284   // bad energy
00285   // ==========
00286   if( !goodenergy ){
00287     return 999.9; // return low resolution
00288   }
00289 
00290   // high resolution
00291   // ===============
00292   if( !fUseBayes 
00293    || !fEnergyPdf
00294    || !fBayesPdf
00295    || ( fCharge==0
00296       && fCosX==0 && fCosY==0 && fCosZ==0
00297        && fEmu==0 && fEshw==0 ) ){
00298     return 0.0;
00299   }
00300 
00301   // bayesan resolution
00302   // ==================
00303   // std::cout << "   using bayesian resolution... " << std::endl;
00304 
00305   // calculate neutrino energy
00306   // =========================
00307   Double_t Enu = fEmu + fEshw;
00308 
00309   // bail out if energy is too high
00310   // ==============================
00311   if( Enu>100.0 ){
00312     return 1.49; // return lowest resolution bin
00313   }
00314 
00315   // loop over neutrino energies
00316   // ===========================
00317   Double_t dEnu = sqrt( 9.0*fEmuErr*fEmuErr + fEshw + 1.0 );
00318   Double_t Emin = Enu - dEnu - 1.0; if( Emin<0.0 ) Emin = 0.0;
00319   Double_t Emax = Enu + dEnu + 2.0; if( Emax>100.0 ) Emax = 100.0;
00320 
00321   Int_t nEmin = (Int_t)(10.0*Emin);    
00322   Int_t nEmax = (Int_t)(10.0*Emax);
00323   Double_t fEtrue = 0.0;
00324 
00325   // std::cout << "   Enu=" << Enu << ": try limits [" << Emin << "," << Emax << "]" << std::endl;
00326 
00327   for( Int_t nE=nEmin; nE<nEmax; nE++ ){
00328     
00329     if( fEnergyWeight[nE]>0.0 ){
00330       fEtrue = fEnergyMedian[nE];
00331       
00332       //
00333       // std::cout << " nE=" << nE << "[" << nEmin << "->" << nEmax << "] Etrue=" << fEtrue << std::endl;
00334       //
00335 
00336       // neutrinos (negative)
00337       // ===================
00338       if( fCharge<=0 ){
00339 
00340         // Quasi-Elastic
00341         this->AddToBayesPdf( fBayesPdf,fNuQE[nE],fEnergyWeight[nE],
00342                              rangecurve,fEtrue,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,true,bayesprior);
00343         
00344         // Resonance
00345         this->AddToBayesPdf( fBayesPdf,fNuRES[nE],fEnergyWeight[nE],
00346                              rangecurve,fEtrue,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,false,bayesprior);
00347 
00348         // Deep Inelastic
00349         this->AddToBayesPdf( fBayesPdf,fNuDIS[nE],fEnergyWeight[nE],
00350                              rangecurve,fEtrue,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,false,bayesprior);
00351       }
00352 
00353       // anti-neutrinos (positive)
00354       // ========================
00355       if( fCharge>=0 ){
00356 
00357         // Quasi-Elastic
00358         this->AddToBayesPdf( fBayesPdf,fNuBarQE[nE],fEnergyWeight[nE],
00359                              rangecurve,fEtrue,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,true,bayesprior);
00360         
00361         // Resonance 
00362         this->AddToBayesPdf( fBayesPdf,fNuBarRES[nE],fEnergyWeight[nE],
00363                              rangecurve,fEtrue,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,false,bayesprior);
00364 
00365         // Deep Inelastic
00366         this->AddToBayesPdf( fBayesPdf,fNuBarDIS[nE],fEnergyWeight[nE],
00367                              rangecurve,fEtrue,fCosX,fCosY,fCosZ,fEmu,fEmuErr,fEshw,false,bayesprior);
00368       }
00369 
00370     } // if( fEnergyWeight[nE]>0.0 ){ ...
00371   } // for( Int_t nEmin=0; nE<nEmax; nE++ ){ ...
00372 
00373   // return rms/mean
00374   // ===============
00375   if( fBayesPdf->GetTotalWeight()>0.0 ){
00376     // std::cout << "   LogLE: Mean=" << fBayesPdf->GetMean() << " Sigma=" << fBayesPdf->GetRMS() << std::endl; 
00377     return fBayesPdf->GetRMS(); // fBayesPdf->GetRMS()/fBayesPdf->GetMean()
00378   }
00379 
00380   // low resolution
00381   // ==============
00382   return 999.9;
00383 }

TH1D * AtNuResolution::GetResolutionPdf ( AtNuResolutionEventType::AtNuResolutionEventType_t  eventtype,
Bool_t  gooddirection,
Bool_t  goodenergy,
Bool_t  rangecurve,
Int_t  fCharge = 0,
Double_t  fCosX = 0.0,
Double_t  fCosY = 0.0,
Double_t  fCosZ = 0.0,
Double_t  fEmu = 0.0,
Double_t  fEmuErr = 0.0,
Double_t  fEshw = 0.0,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior = AtNuResolutionPrior::kFlat,
Bool_t  renormalise = 0 
)

Definition at line 154 of file AtNuResolution.cxx.

References fBayesPdf, fHistNumber, GetResolution(), AtNuResolutionEventType::kCV, AtNuResolutionEventType::kNUE, AtNuResolutionEventType::kUPMU, and AtNuResBayesPdf::MakeTH1D().

00155 {
00156   TH1D* hResolution = 0;
00157 
00158   GetResolution(eventtype, 
00159                 gooddirection, goodenergy, rangecurve, 
00160                 fCharge, fCosX, fCosY, fCosZ, fEmu, fEmuErr, fEshw,
00161                 bayesprior);
00162 
00163   switch( eventtype ){
00164     case AtNuResolutionEventType::kNUE: 
00165       break;
00166     case AtNuResolutionEventType::kCV:
00167       if( fBayesPdf ){
00168         TString histname("bayespdf");
00169         if( fHistNumber<10 ) histname.Append("000");
00170         else if( fHistNumber<100 ) histname.Append("00");
00171         else if( fHistNumber<1000 ) histname.Append("0");
00172         histname += fHistNumber; fHistNumber++;
00173         hResolution = (TH1D*)fBayesPdf->MakeTH1D(histname.Data(),renormalise); 
00174       }
00175       break;
00176     case AtNuResolutionEventType::kUPMU:
00177       break;
00178     default:
00179       break;
00180   }
00181 
00182   return hResolution;
00183 }

AtNuResolutionType::AtNuResolutionType_t AtNuResolution::GetResolutionType ( AtNuResolutionEventType::AtNuResolutionEventType_t  eventtype,
Double_t  resolution 
)

Definition at line 195 of file AtNuResolution.cxx.

References AtNuResolutionEventType::kCV, AtNuResolutionType::kHiRes1, AtNuResolutionType::kHiRes2, AtNuResolutionType::kHiRes3, AtNuResolutionType::kHiRes4, AtNuResolutionType::kLoRes, AtNuResolutionType::kNone, AtNuResolutionEventType::kNUE, and AtNuResolutionEventType::kUPMU.

00196 {
00197   switch( eventtype ){
00198     case AtNuResolutionEventType::kNUE: 
00199       return AtNuResolutionType::kLoRes;
00200     case AtNuResolutionEventType::kCV:
00201       if( resolution<0.25 ) return AtNuResolutionType::kHiRes4;
00202       else if( resolution<0.50 ) return AtNuResolutionType::kHiRes3;
00203       else if( resolution<0.75 ) return AtNuResolutionType::kHiRes2;
00204       else if( resolution<1.50 ) return AtNuResolutionType::kHiRes1;
00205       else return AtNuResolutionType::kLoRes;
00206     case AtNuResolutionEventType::kUPMU:
00207       if( resolution<0.25 ) return AtNuResolutionType::kHiRes4;
00208       else if( resolution<0.50 ) return AtNuResolutionType::kHiRes3;
00209       else if( resolution<0.75 ) return AtNuResolutionType::kHiRes2;
00210       else if( resolution<1.50 ) return AtNuResolutionType::kHiRes1;
00211       else return AtNuResolutionType::kLoRes;
00212       return AtNuResolutionType::kLoRes;
00213     default:
00214       return AtNuResolutionType::kNone;
00215   }
00216 
00217   return AtNuResolutionType::kNone;
00218 }

AtNuResolutionType::AtNuResolutionType_t AtNuResolution::GetResolutionType ( AtNuResolutionEventType::AtNuResolutionEventType_t  eventtype,
Bool_t  gooddirection,
Bool_t  goodenergy,
Bool_t  rangecurve,
Int_t  fCharge = 0,
Double_t  fCosX = 0.0,
Double_t  fCosY = 0.0,
Double_t  fCosZ = 0.0,
Double_t  fEmu = 0.0,
Double_t  fEmuErr = 0.0,
Double_t  fEshw = 0.0,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior = AtNuResolutionPrior::kFlat 
)

Definition at line 185 of file AtNuResolution.cxx.

References GetResolution().

Referenced by AtNuResNtuple::Run().

00186 {
00187   Double_t resolution = GetResolution(eventtype, 
00188                                       gooddirection, goodenergy, rangecurve,
00189                                       fCharge, fCosX, fCosY, fCosZ, fEmu, fEmuErr, fEshw,
00190                                       bayesprior);
00191   
00192   return GetResolutionType(eventtype,resolution);
00193 }  

Double_t AtNuResolution::GetResolutionUPMU ( Bool_t  gooddirection,
Bool_t  goodenergy,
Bool_t  rangecurve,
Int_t  fCharge = 0,
Double_t  fCosX = 0.0,
Double_t  fCosY = 0.0,
Double_t  fCosZ = 0.0,
Double_t  fEmu = 0.0,
Double_t  fEmuErr = 0.0,
Double_t  fEshw = 0.0,
AtNuResolutionPrior::AtNuResolutionPrior_t  bayesprior = AtNuResolutionPrior::kFlat 
) [private]

Definition at line 220 of file AtNuResolution.cxx.

References fBayesPdf, fEnergyPdf, fUseBayes, and AtNuResBayesPdf::Reset().

Referenced by GetResolution().

00221 {
00222   // Upward Muon Events
00223   // ==================
00224   // std::cout << " *** AtNuResolution::GetResolutionUPMU(...) *** " << std::endl;
00225   // std::cout << "   gooddirection=" << gooddirection << ", goodenergy=" << goodenergy << ", rangecurve=" << rangecurve << std::endl << "   charge=" << fCharge << ", cosX=" << fCosX << ", cosY=" << fCosY << ", cosZ=" << fCosZ << ", emu=" << fEmu << ", emuerr=" << fEmuErr << ", eshw=" << fEshw << std::endl;
00226 
00227   // Reset Bayes Pdf
00228   // ===============
00229   if( fBayesPdf ) fBayesPdf->Reset();
00230 
00231   // bad direction
00232   // =============
00233   if( !gooddirection ){
00234     return 999.9; // return low resolution
00235   }
00236 
00237   // bad energy
00238   // ==========
00239   if( !goodenergy ){
00240     return 0.45; // return high momentum
00241   }
00242 
00243   // high resolution
00244   // ===============
00245   if( !fUseBayes 
00246    || !fEnergyPdf
00247    || !fBayesPdf
00248    || ( fCharge==0
00249       && fCosX==0 && fCosY==0 && fCosZ==0
00250        && fEmu==0 && fEshw==0 ) ){
00251     return 0.0;
00252   }
00253 
00254   // separate by muon momentum 
00255   // =========================
00256   // place momentum cut at 10 GeV
00257   if( fEmu<10.0 ) return 0.15; 
00258   else return 0.30;
00259   
00260   // low resolution
00261   // ==============
00262   return 999.9;
00263 }

void AtNuResolution::LoadHistograms ( const char *  filename,
Bool_t  flag = 1 
)

Definition at line 597 of file AtNuResolution.cxx.

References LoadHistogramsFromCambridge().

Referenced by AtNuResNtuple::LoadHistograms().

00598 {
00599   std::cout << " *** AtNuResolution::LoadHistograms(...) *** " << std::endl;  if( flag ){
00600     std::cout << "   Loading histograms from Cambridge... " << std::endl;
00601     this->LoadHistogramsFromCambridge(filename);
00602   }
00603 
00604   return;
00605 }

void AtNuResolution::LoadHistogramsFromCambridge ( const char *  filename  )  [private]

Definition at line 607 of file AtNuResolution.cxx.

References fEnergyMedian, fEnergyWeight, fNeugenPdf, fNuBarDIS, fNuBarQE, fNuBarRES, fNuDIS, fNuQE, fNuRES, fUseBayes, and n.

Referenced by LoadHistograms().

00608 {
00609   std::cout << " *** AtNuResolution::LoadHistogramsFromCambridge(...) *** " << std::endl
00610             << "   Loading histograms from: " << filename << std::endl;
00611 
00612   TDirectory* tmpd = gDirectory;
00613   TFile* myFile = TFile::Open(filename,"read");
00614 
00615   if( myFile ){
00616 
00617     TString histname = "";
00618     TString appendname = "";
00619 
00620     TH1D* h1 = 0;
00621     TH2D* h2 = 0;
00622 
00623     Int_t bin = 0;
00624     Double_t weight = 0.0;
00625     Double_t energyweight = 0.0;
00626 
00627     TH1D* hNuQEfrac = (TH1D*)(myFile->Get("NuQEfrac"));
00628     TH1D* hNuResfrac = (TH1D*)(myFile->Get("NuResfrac"));
00629     TH1D* hNuDisfrac = (TH1D*)(myFile->Get("NuDisfrac"));
00630 
00631     TH1D* hNuBarQEfrac = (TH1D*)(myFile->Get("NuBarQEfrac"));
00632     TH1D* hNuBarResfrac = (TH1D*)(myFile->Get("NuBarResfrac"));
00633     TH1D* hNuBarDisfrac = (TH1D*)(myFile->Get("NuBarDisfrac"));
00634 
00635     for( Int_t n=0; n<278; n++ ){
00636 
00637       appendname = "";
00638       Int_t lowEnergy  = 0;
00639       Int_t highEnergy = 0;
00640 
00641       if( n<98 ){
00642         lowEnergy  += 200+100*(n);
00643         highEnergy += 200+100*(n+1);
00644       }
00645       else{
00646         lowEnergy  += 10000+500*(n-98);
00647         highEnergy += 10000+500*(n-98+1);
00648       }
00649 
00650       appendname += lowEnergy;
00651       appendname.Append("MeV");
00652       appendname += highEnergy;
00653       appendname.Append("MeV");
00654 
00655       bin = (Int_t)(lowEnergy/100.0);
00656       energyweight = (highEnergy-lowEnergy)/500.0;
00657 
00658       // neutrino QE
00659       // ===========
00660       histname = "NuQEy";
00661       histname.Append(appendname.Data());
00662       weight = hNuQEfrac->GetBinContent(n+1);
00663 
00664       h1 = (TH1D*)(myFile->Get(histname.Data()));
00665       fNeugenPdf = new AtNuResNeugenPdf(h1,+14,weight);
00666       fNuQE[bin] = fNeugenPdf;
00667 
00668       // neutrino RES
00669       // ============
00670       histname = "NuResonanceYvsW2";
00671       histname.Append(appendname.Data());
00672       weight = hNuResfrac->GetBinContent(n+1);
00673 
00674       h2 = (TH2D*)(myFile->Get(histname.Data()));
00675       fNeugenPdf = new AtNuResNeugenPdf(h2,+14,weight);
00676       fNuRES[bin] = fNeugenPdf;
00677 
00678       // neutrino DIS
00679       // ============
00680       histname = "NuDisYvsW2";
00681       histname.Append(appendname.Data());
00682       weight = hNuDisfrac->GetBinContent(n+1);
00683 
00684       h2 = (TH2D*)(myFile->Get(histname.Data()));
00685       fNeugenPdf = new AtNuResNeugenPdf(h2,+14,weight);
00686       fNuDIS[bin] = fNeugenPdf; 
00687 
00688       // anti-neutrino QE
00689       // ================
00690       histname = "NuBarQEy";
00691       histname.Append(appendname.Data());
00692       weight = hNuBarQEfrac->GetBinContent(n+1);
00693 
00694       h1 = (TH1D*)(myFile->Get(histname.Data()));
00695       fNeugenPdf = new AtNuResNeugenPdf(h1,-14,weight);
00696       fNuBarQE[bin] = fNeugenPdf;
00697 
00698       // anti-neutrino RES
00699       // =================
00700       histname = "NuBarResonanceYvsW2";
00701       histname.Append(appendname.Data());
00702       weight = hNuBarResfrac->GetBinContent(n+1);
00703 
00704       h2 = (TH2D*)(myFile->Get(histname.Data()));
00705       fNeugenPdf = new AtNuResNeugenPdf(h2,-14,weight);
00706       fNuBarRES[bin] = fNeugenPdf;
00707 
00708       // anti-neutrino DIS
00709       // =================
00710       histname = "NuBarDisYvsW2";
00711       histname.Append(appendname.Data());
00712       weight = hNuBarDisfrac->GetBinContent(n+1);
00713 
00714       h2 = (TH2D*)(myFile->Get(histname.Data()));
00715       fNeugenPdf = new AtNuResNeugenPdf(h2,-14,weight);
00716       fNuBarDIS[bin] = fNeugenPdf;
00717 
00718       // energy weighting
00719       // ================
00720       fEnergyWeight[bin] = energyweight;
00721       fEnergyMedian[bin] = 0.5*0.001*(lowEnergy+highEnergy);
00722     }
00723 
00724     fUseBayes = 1;
00725   }
00726 
00727   myFile->Close();
00728   gDirectory = tmpd;
00729 
00730   return;
00731 }

void AtNuResolution::SetEnergyBinning ( Int_t  bins,
Double_t  min,
Double_t  max 
)

Definition at line 587 of file AtNuResolution.cxx.

References fBayesPdf.

Referenced by AtNuResNtuple::SetEnergyBinning().

00588 {
00589   std::cout << " *** AtNuResolution::SetEnergyBinning(...) *** " << std::endl;
00590   std::cout << "   binning = [" << bins << "," << min << "->" << max << "]" << std::endl; 
00591 
00592   if( !fBayesPdf ){
00593     fBayesPdf = new AtNuResBayesPdf(bins,min,max);
00594   }
00595 }


Member Data Documentation

Double_t* AtNuResolution::fEnergyMedian [private]

Definition at line 224 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Double_t* AtNuResolution::fEnergyWeight [private]

Definition at line 223 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Int_t AtNuResolution::fHistNumber [private]

Definition at line 227 of file AtNuResolution.h.

Referenced by GetResolutionPdf().

Definition at line 213 of file AtNuResolution.h.

Referenced by LoadHistogramsFromCambridge().

Definition at line 221 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Definition at line 219 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Definition at line 220 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Definition at line 217 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Definition at line 215 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Definition at line 216 of file AtNuResolution.h.

Referenced by GetResolutionCV(), LoadHistogramsFromCambridge(), and ~AtNuResolution().

Bool_t AtNuResolution::fUseBayes [private]

The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1