loon: NuZBeamReweight.cxx Source File

00001 
00002 
00003 // Coded by Jeff Hartnell July/2006       
00004 //                                                                    
00005 // Contact: j.j.hartnell@rl.ac.uk
00007 
00008 #include "TClonesArray.h"
00009 #include "TDirectory.h"
00010 
00011 #include "Conventions/BeamType.h"
00012 #include "MessageService/MsgService.h"
00013 #include "MCNtuple/NtpMCTruth.h"
00014 #include "CandNtupleSR/NtpSREvent.h"
00015 #include "TruthHelperNtuple/NtpTHEvent.h"
00016 #include "StandardNtuple/NtpStRecord.h"
00017 #include "Conventions/ReleaseType.h"
00018 
00019 #include "MCReweight/MCReweight.h"
00020 #include "MCReweight/NeugenWeightCalculator.h"
00021 #include "MCReweight/ReweightHelpers.h"
00022 #include "MCReweight/SKZPWeightCalculator.h"
00023 #include "MCReweight/Zbeam.h"
00024 #include "MCReweight/Zfluk.h"
00025 
00026 #include "NtupleUtils/NuEvent.h"
00027 #include "NtupleUtils/NuZBeamReweight.h"
00028 
00029 CVSID("$Id: NuZBeamReweight.cxx,v 1.25 2008/02/20 14:34:23 hartnell Exp $");
00030 
00031 //......................................................................
00032 
00033 NuZBeamReweight::NuZBeamReweight()
00034 {
00035   MSG("NuZBeamReweight",Msg::kDebug)
00036     <<"Running NuZBeamReweight Constructor..."<<endl;
00037 
00038 
00039   MSG("NuZBeamReweight",Msg::kDebug)
00040     <<"Finished NuZBeamReweight Constructor"<<endl;
00041 }
00042 
00043 //......................................................................
00044 
00045 NuZBeamReweight::~NuZBeamReweight()
00046 {
00047   MSG("NuZBeamReweight",Msg::kDebug)
00048     <<"Running NuZBeamReweight Destructor..."<<endl;
00049   
00050 
00051   MSG("NuZBeamReweight",Msg::kDebug)
00052     <<"Finished NuZBeamReweight Destructor"<<endl;
00053 }
00054 
00055 //......................................................................
00056 
00057 void NuZBeamReweight::ExtractZBeamReweight(NuEvent& nu) const
00058 {
00059   if (nu.simFlag!=SimFlag::kMC || 
00060       nu.reweightVersion==SKZPWeightCalculator::kUnknown) {
00061     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00062       <<"Not doing ExtractZBeamReweight, reweightVersion="
00063       <<nu.reweightVersion<<", simFlag="<<nu.simFlag<<endl;
00064     
00065     //set the weights to 1 that would otherwise be set below
00066     nu.trkEnWeight=1;//no reweighting
00067     nu.shwEnWeight=1;//no reweighting
00068     nu.beamWeight=1;//no reweighting
00069     nu.fluxErrHadProdAfterTune=-1;
00070     nu.fluxErrTotalErrorPreTune=-1;
00071     nu.fluxErrTotalErrorAfterTune=-1;
00072     nu.detectorWeightNMB=1;//no reweighting
00073     nu.detectorWeightNM=1;//no reweighting
00074     
00075     nu.trkEnWeightRunI=1;//no reweighting
00076     nu.shwEnWeightRunI=1;//no reweighting
00077     nu.beamWeightRunI=1;//no reweighting
00078     nu.fluxErrHadProdAfterTuneRunI=-1;
00079     nu.fluxErrTotalErrorPreTuneRunI=-1;
00080     nu.fluxErrTotalErrorAfterTuneRunI=-1;
00081     nu.detectorWeightNMBRunI=1;//no reweighting
00082     nu.detectorWeightNMRunI=1;//no reweighting
00083     
00084     nu.trkEnWeightRunII=1;//no reweighting
00085     nu.shwEnWeightRunII=1;//no reweighting
00086     nu.beamWeightRunII=1;//no reweighting
00087     nu.fluxErrHadProdAfterTuneRunII=-1;
00088     nu.fluxErrTotalErrorPreTuneRunII=-1;
00089     nu.fluxErrTotalErrorAfterTuneRunII=-1;
00090     nu.detectorWeightNMBRunII=1;//no reweighting
00091     nu.detectorWeightNMRunII=1;//no reweighting
00092 
00093     return;
00094   }
00095 
00096   MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00097     <<"Running ExtractZBeamReweight..."<<endl;
00098 
00099   //convert the std BeamType into Zarko's convention
00100   //was 2=le010z185i (Zbeam.h)
00101   const Int_t Ibeam=BeamType::ToZarko
00102     (static_cast<BeamType::BeamType_t>(nu.beamType));
00103     
00104   //get a reference
00105   SKZPWeightCalculator& wc=this->GetSKZPWeightCalculator(nu);
00106     
00107   //copied from SKZPWeightCalculator
00108   Int_t ccnc=nu.iaction;
00109   Double_t true_enu=nu.neuEnMC;//p4neu[3]
00110   Double_t true_eshw=nu.shwEnMC;//p4shw[3]
00111   Int_t inu=nu.inu; 
00112   Float_t nsigma=1;
00113   //FUTURE VERSIONS may use this:
00114   //parenttype = fi->ptype;
00115   //parentpz = 1.*fi->pppz;
00116   //parentpy = 1.*fi->ppdydz*parentpz;
00117   //parentpx = 1.*fi->ppdxdz*parentpz;
00118   Int_t parenttype = nu.tptype;
00119   Double_t parentpz = nu.tpz;
00120   Double_t parentpy = nu.tpy;
00121   Double_t parentpx = nu.tpx;
00122   Double_t pt = sqrt(parentpy*parentpy+parentpx*parentpx);
00123 
00124   //fill MCEventInfo
00125   MCEventInfo ei;
00126   ei.UseStoredXSec(true);
00127   this->MCEventInfoFilla(&ei,nu);
00128 
00130   //calc the weights for runI
00132   SKZPWeightCalculator::RunPeriod_t rp=SKZPWeightCalculator::kRunI;
00133 
00134   //create variables that are passed by reference to function to fill
00135   Double_t new_reco_emu1=0;
00136   Double_t new_reco_eshw1=0;
00137   Double_t beamweight1=0;
00138   Double_t detweightNMB1=0;
00139   Double_t detweightNM1=0;
00140     
00141   //get the skzp weights (beam and detector)
00142   wc.SetSampleSelection("NuMuBar");
00143   wc.GetWeight(nu.detector,Ibeam,
00144                parenttype,pt,parentpz,ccnc,true_enu,inu,
00145                nu.trkEn,nu.shwEn,new_reco_emu1,new_reco_eshw1,
00146                beamweight1,detweightNMB1,rp,true_eshw,&ei);
00147 
00148   MAXMSG("NuZBeamReweight",Msg::kDebug,200)
00149     <<endl<<"iaction="<<nu.iaction
00150     <<", nu.trkEn="<<nu.trkEn<<", nu.shwEn="<<nu.shwEn<<endl;
00151 
00152   MAXMSG("NuZBeamReweight",Msg::kDebug,200)
00153     <<"NMB: newEmu1="<<new_reco_emu1
00154     <<", newEShw1="<<new_reco_eshw1
00155     <<", detNMB="<<detweightNMB1
00156     <<", beamW="<<beamweight1<<endl;
00157 
00158   //get the skzp weights (beam and detector)
00159   wc.SetSampleSelection("NuMu");
00160   wc.GetWeight(nu.detector,Ibeam,
00161                parenttype,pt,parentpz,ccnc,true_enu,inu,
00162                nu.trkEn,nu.shwEn,new_reco_emu1,new_reco_eshw1,
00163                beamweight1,detweightNM1,rp,true_eshw,&ei);
00164   
00165   MAXMSG("NuZBeamReweight",Msg::kDebug,200)
00166     <<"NM:  newEmu1="<<new_reco_emu1
00167     <<", newEshw1="<<new_reco_eshw1
00168     <<",  detNM="<<detweightNM1
00169     <<", beamW="<<beamweight1<<endl;
00170 
00171   //new for kDetXs
00172   //double GetWeight(int det, int Ibeam,
00173   //       int tptype, double pt, double pz,    
00175   //       double reco_emu, double reco_eshw,
00176   //       double &new_reco_emu, double &new_reco_eshw,
00177   //       double &beamweight, double &detweight,
00178   //       RunPeriod_t rp=kRunI,double true_eshw=-1,MCEventInfo *ei=0);
00179 
00180   Float_t trkEnWeight1=1;
00181   Float_t shwEnWeight1=1;
00182   if (nu.trkEn) trkEnWeight1=new_reco_emu1/nu.trkEn;
00183   if (nu.shwEn) shwEnWeight1=new_reco_eshw1/nu.shwEn;
00184 
00185   //get the flux errors
00186   Float_t fluxErrHadProdAfterTune1=wc.GetFluxError
00187     (nu.detector,Ibeam,inu,true_enu,
00188      SKZPWeightCalculator::kHadProdAfterTune,nsigma);
00189   Float_t fluxErrTotalErrorPreTune1=wc.GetFluxError
00190     (nu.detector,Ibeam,inu,true_enu,
00191      SKZPWeightCalculator::kTotalErrorPreTune,nsigma);
00192   Float_t fluxErrTotalErrorAfterTune1=wc.GetFluxError
00193     (nu.detector,Ibeam,inu,true_enu,
00194      SKZPWeightCalculator::kTotalErrorAfterTune,nsigma);
00195 
00197   //calc the weights for runII
00199   rp=SKZPWeightCalculator::kRunII;
00200     
00201   //create variables that are passed by reference to function to fill
00202   Double_t new_reco_emu2=0;
00203   Double_t new_reco_eshw2=0;
00204   Double_t beamweight2=0;
00205   Double_t detweightNMB2=0;
00206   Double_t detweightNM2=0;
00207     
00208   //get the skzp weights (beam and detector)
00209   wc.SetSampleSelection("NuMuBar");
00210   wc.GetWeight(nu.detector,Ibeam,
00211                parenttype,pt,parentpz,ccnc,true_enu,inu,
00212                nu.trkEn,nu.shwEn,new_reco_emu2,new_reco_eshw2,
00213                beamweight2,detweightNMB2,rp,true_eshw,&ei);
00214 
00215   //get the skzp weights (beam and detector)
00216   wc.SetSampleSelection("NuMu");
00217   wc.GetWeight(nu.detector,Ibeam,
00218                parenttype,pt,parentpz,ccnc,true_enu,inu,
00219                nu.trkEn,nu.shwEn,new_reco_emu2,new_reco_eshw2,
00220                beamweight2,detweightNM2,rp,true_eshw,&ei);
00221 
00222   Float_t trkEnWeight2=1;
00223   Float_t shwEnWeight2=1;
00224   if (nu.trkEn) trkEnWeight2=new_reco_emu2/nu.trkEn;
00225   if (nu.shwEn) shwEnWeight2=new_reco_eshw2/nu.shwEn;
00226     
00227   //get the flux errors
00228   Float_t fluxErrHadProdAfterTune2=wc.GetFluxError
00229     (nu.detector,Ibeam,inu,true_enu,
00230      SKZPWeightCalculator::kHadProdAfterTune,nsigma);
00231   Float_t fluxErrTotalErrorPreTune2=wc.GetFluxError
00232     (nu.detector,Ibeam,inu,true_enu,
00233      SKZPWeightCalculator::kTotalErrorPreTune,nsigma);
00234   Float_t fluxErrTotalErrorAfterTune2=wc.GetFluxError
00235     (nu.detector,Ibeam,inu,true_enu,
00236      SKZPWeightCalculator::kTotalErrorAfterTune,nsigma);
00237 
00239   //calc the weights for the weighted average
00241   static Bool_t firstTime=true;
00242   if (firstTime) {
00243     firstTime=false;
00244     //define and set the numbers of POTs in each run period
00245     //numbers from Andy Blakes minos-doc-3727
00246     Float_t potRunI=1.269e20;
00247     //runIIa=1.229
00248     //runIIb=0.721
00249     Float_t potRunII=1.950e20;
00250     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00251       <<"ExtractZBeamReweight: SetRunFrac w/ potRunI="<<potRunI
00252       <<", potRunII="<<potRunII<<endl;
00253     wc.SetRunFrac(Ibeam,SKZPWeightCalculator::kRunI,potRunI);
00254     wc.SetRunFrac(Ibeam,SKZPWeightCalculator::kRunII,potRunII);
00255   }
00256     
00257   //create variables that are passed by reference to function to fill
00258   Double_t new_reco_emuAv=0;
00259   Double_t new_reco_eshwAv=0;
00260   Double_t beamweightAv=0;
00261   Double_t detweightNMBAv=0;
00262   Double_t detweightNMAv=0;
00263     
00264   //get the skzp weights (beam and detector)
00265   wc.SetSampleSelection("NuMuBar");
00266   wc.GetRFWWeight(nu.detector,Ibeam,
00267                   parenttype,pt,parentpz,ccnc,true_enu,inu,
00268                   nu.trkEn,nu.shwEn,new_reco_emuAv,new_reco_eshwAv,
00269                   beamweightAv,detweightNMBAv);
00270 
00271   //get the skzp weights (beam and detector)
00272   wc.SetSampleSelection("NuMu");
00273   wc.GetRFWWeight(nu.detector,Ibeam,
00274                   parenttype,pt,parentpz,ccnc,true_enu,inu,
00275                   nu.trkEn,nu.shwEn,new_reco_emuAv,new_reco_eshwAv,
00276                   beamweightAv,detweightNMAv);
00277 
00278   Float_t trkEnWeightAv=1;
00279   Float_t shwEnWeightAv=1;
00280   if (nu.trkEn) trkEnWeightAv=new_reco_emuAv/nu.trkEn;
00281   if (nu.shwEn) shwEnWeightAv=new_reco_eshwAv/nu.shwEn;
00282     
00283   //get the flux errors
00284   Float_t fluxErrHadProdAfterTune=wc.GetFluxError
00285     (nu.detector,Ibeam,inu,true_enu,
00286      SKZPWeightCalculator::kHadProdAfterTune,nsigma);
00287   Float_t fluxErrTotalErrorPreTune=wc.GetFluxError
00288     (nu.detector,Ibeam,inu,true_enu,
00289      SKZPWeightCalculator::kTotalErrorPreTune,nsigma);
00290   Float_t fluxErrTotalErrorAfterTune=wc.GetFluxError
00291     (nu.detector,Ibeam,inu,true_enu,
00292      SKZPWeightCalculator::kTotalErrorAfterTune,nsigma);
00293 
00294   MAXMSG("NuZBeamReweight",Msg::kDebug,300)
00295     <<"Wgts RunI:  b="<<beamweight1
00296     <<", dNMB="<<detweightNMB1<<", dNM="<<detweightNM1
00297     <<", trk="<<trkEnWeight1<<", shw="<<shwEnWeight1<<endl
00298     <<"Wgts RunII: b="<<beamweight2
00299     <<", dNMB="<<detweightNMB2<<", dNM="<<detweightNM2
00300     <<", trk="<<trkEnWeight2<<", shw="<<shwEnWeight2<<endl
00301     <<"Wgts Av:    b="<<beamweightAv
00302     <<", dNMB="<<detweightNMBAv<<", dNM="<<detweightNMAv
00303     <<", trk="<<trkEnWeightAv<<", shw="<<shwEnWeightAv<<endl<<endl;
00304   
00305   MAXMSG("NuZBeamReweight",Msg::kDebug,200)
00306     <<"energyMC="<<nu.neuEnMC
00307     <<", iaction="<<nu.iaction
00308     <<", inu="<<nu.inu
00309     <<endl
00310     <<"fluxErrs:"<<" HadProdAfterTune1   ="
00311     <<fluxErrHadProdAfterTune1<<endl
00312     <<"         "<<" TotalErrorPreTune1  ="
00313     <<fluxErrTotalErrorPreTune1<<endl
00314     <<"         "<<" TotalErrorAfterTune1="
00315     <<fluxErrTotalErrorAfterTune1<<endl
00316     <<"fluxErrs:"<<" HadProdAfterTune2   ="
00317     <<fluxErrHadProdAfterTune2<<endl
00318     <<"         "<<" TotalErrorPreTune2  ="
00319     <<fluxErrTotalErrorPreTune2<<endl
00320     <<"         "<<" TotalErrorAfterTune2="
00321     <<fluxErrTotalErrorAfterTune2<<endl
00322     <<"fluxErrs:"<<" HadProdAfterTune   ="
00323     <<fluxErrHadProdAfterTune<<endl
00324     <<"         "<<" TotalErrorPreTune  ="
00325     <<fluxErrTotalErrorPreTune<<endl
00326     <<"         "<<" TotalErrorAfterTune="
00327     <<fluxErrTotalErrorAfterTune<<endl;
00328  
00330   //store the weights
00332   nu.trkEnWeight=trkEnWeightAv;
00333   nu.shwEnWeight=shwEnWeightAv;
00334   nu.beamWeight=beamweightAv;
00335   nu.fluxErrHadProdAfterTune=fluxErrHadProdAfterTune;
00336   nu.fluxErrTotalErrorPreTune=fluxErrTotalErrorPreTune;
00337   nu.fluxErrTotalErrorAfterTune=fluxErrTotalErrorAfterTune;
00338   nu.detectorWeightNMB=detweightNMBAv;
00339   nu.detectorWeightNM=detweightNMAv;
00340     
00341   nu.trkEnWeightRunI=trkEnWeight1;
00342   nu.shwEnWeightRunI=shwEnWeight1;
00343   nu.beamWeightRunI=beamweight1;
00344   nu.fluxErrHadProdAfterTuneRunI=fluxErrHadProdAfterTune1;
00345   nu.fluxErrTotalErrorPreTuneRunI=fluxErrTotalErrorPreTune1;
00346   nu.fluxErrTotalErrorAfterTuneRunI=fluxErrTotalErrorAfterTune1;
00347   nu.detectorWeightNMBRunI=detweightNMB1;
00348   nu.detectorWeightNMRunI=detweightNM1;
00349     
00350   nu.trkEnWeightRunII=trkEnWeight2;
00351   nu.shwEnWeightRunII=shwEnWeight2;
00352   nu.beamWeightRunII=beamweight2;
00353   nu.fluxErrHadProdAfterTuneRunII=fluxErrHadProdAfterTune2;
00354   nu.fluxErrTotalErrorPreTuneRunII=fluxErrTotalErrorPreTune2;
00355   nu.fluxErrTotalErrorAfterTuneRunII=fluxErrTotalErrorAfterTune2;
00356   nu.detectorWeightNMBRunII=detweightNMB2;
00357   nu.detectorWeightNMRunII=detweightNM2;
00358 }
00359 
00360 //......................................................................
00361 
00362 void NuZBeamReweight::ExtractZBeamReweight(const NtpStRecord& /*ntp*/,
00363                                            const NtpSREvent& /*evt*/,
00364                                            NuEvent& nu) const
00365 { 
00366   //get the weights
00367   this->ExtractZBeamReweight(nu);
00368 }
00369 
00370 //......................................................................
00371 
00372 void NuZBeamReweight::ExtractZBeamReweight(const NtpStRecord& ntp,
00373                                            NuEvent& nu) const
00374 { 
00375   MAXMSG("NuZBeamReweight",Msg::kError,2000)
00376     <<"NuZBeamReweight::ExtractZBeamReweight(const NtpStRecord& ntp,NuEvent& nu) const"
00377     <<" is deprecated!"<<endl;
00378   
00379   MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00380     <<"Running ExtractZBeamReweight..."<<endl;
00381 
00382   //convert the std BeamType into Zarko's convention
00383   //was 2=le010z185i (Zbeam.h)
00384   const Int_t Ibeam=BeamType::ToZarko
00385     (static_cast<BeamType::BeamType_t>(nu.beamType));
00386     
00387   //get a reference
00388   SKZPWeightCalculator& wc=this->GetSKZPWeightCalculator(nu);
00389     
00390   //copied from SKZPWeightCalculator
00391   Int_t ccnc=nu.iaction;
00392   Double_t true_enu=nu.neuEnMC;//p4neu[3]
00393   Double_t true_eshw=nu.shwEnMC;//p4shw[3]
00394   Int_t inu=nu.inu; 
00395   //Float_t nsigma=1;
00396   //FUTURE VERSIONS may use this:
00397   //parenttype = fi->ptype;
00398   //parentpz = 1.*fi->pppz;
00399   //parentpy = 1.*fi->ppdydz*parentpz;
00400   //parentpx = 1.*fi->ppdxdz*parentpz;
00401   Int_t parenttype = nu.tptype;
00402   Double_t parentpz = nu.tpz;
00403   Double_t parentpy = nu.tpy;
00404   Double_t parentpx = nu.tpx;
00405   Double_t pt = sqrt(parentpy*parentpy+parentpx*parentpx);
00406 
00407   //fill MCEventInfo
00408   //MCEventInfo ei;
00409   //ei.UseStoredXSec(true);
00410   //this->MCEventInfoFilla(&ei,nu);
00411 
00412   //fill MCEventInfo to do MODBYRS3;
00413   MCEventInfo ei;
00414   ei.UseStoredXSec(true);
00415   NtpStRecord* st=const_cast<NtpStRecord*>(&ntp);//nasty!
00416   ReweightHelpers::MCEventInfoFilla(&ei,st,nu.mc);
00417 
00418   MAXMSG("NuZBeamReweight",Msg::kInfo,100)
00419     <<endl<<endl
00420     <<", ei.nuE="<<ei.nuE
00421     <<", ei.nuPx="<<ei.nuPx
00422     <<", ei.nuPy="<<ei.nuPy
00423     <<", ei.nuPz="<<ei.nuPz
00424     <<endl
00425     <<", ei.tarE="<<ei.tarE
00426     <<", ei.tarPx="<<ei.tarPx
00427     <<", ei.tarPy="<<ei.tarPy
00428     <<", ei.tarPz="<<ei.tarPz
00429     <<endl
00430     <<", ei.y="<<ei.y<<endl
00431     <<", ei.x="<<ei.x<<endl
00432     <<", ei.q2="<<ei.q2<<endl
00433     <<", ei.w2="<<ei.w2<<endl
00434     <<endl
00435     <<", ei.iaction="<<ei.iaction<<endl
00436     <<", ei.inu="<<ei.inu<<endl
00437     <<", ei.iresonance="<<ei.iresonance<<endl
00438     <<", ei.initial_state="<<ei.initial_state<<endl
00439     <<", ei.nucleus="<<ei.nucleus<<endl
00440     <<", ei.had_fs="<<ei.had_fs<<endl
00441     <<endl;
00442 
00444   //calc the weights for runI
00446   SKZPWeightCalculator::RunPeriod_t rp=SKZPWeightCalculator::kRunI;
00447 
00448   //create variables that are passed by reference to function to fill
00449   Double_t new_reco_emu1=0;
00450   Double_t new_reco_eshw1=0;
00451   Double_t beamweight1=0;
00452   Double_t detweightNMB1=0;
00453   Double_t detweightNM1=0;
00454     
00455   //get the skzp weights (beam and detector)
00456   wc.SetSampleSelection("NuMuBar");
00457   wc.GetWeight(nu.detector,Ibeam,
00458                parenttype,pt,parentpz,ccnc,true_enu,inu,
00459                nu.trkEn,nu.shwEn,new_reco_emu1,new_reco_eshw1,
00460                beamweight1,detweightNMB1,rp,true_eshw,&ei);
00461 
00462   //get the skzp weights (beam and detector)
00463   wc.SetSampleSelection("NuMu");
00464   wc.GetWeight(nu.detector,Ibeam,
00465                parenttype,pt,parentpz,ccnc,true_enu,inu,
00466                nu.trkEn,nu.shwEn,new_reco_emu1,new_reco_eshw1,
00467                beamweight1,detweightNM1,rp,true_eshw,&ei);
00468   
00469   //new for kDetXs
00470   //double GetWeight(int det, int Ibeam,
00471   //       int tptype, double pt, double pz,    
00473   //       double reco_emu, double reco_eshw,
00474   //       double &new_reco_emu, double &new_reco_eshw,
00475   //       double &beamweight, double &detweight,
00476   //       RunPeriod_t rp=kRunI,double true_eshw=-1,MCEventInfo *ei=0);
00477 
00478   Float_t trkEnWeight1=1;
00479   Float_t shwEnWeight1=1;
00480   if (nu.trkEn) trkEnWeight1=new_reco_emu1/nu.trkEn;
00481   if (nu.shwEn) shwEnWeight1=new_reco_eshw1/nu.shwEn;
00482 }
00483 
00484 //......................................................................
00485 
00486 SKZPWeightCalculator& NuZBeamReweight::
00487 GetSKZPWeightCalculator(const NuEvent& nu) const
00488 {
00489   //make a weight calculator, defining the skzp version to use
00490   static SKZPWeightCalculator* pwc=0;
00491   
00492   //create the weight calculator on the first call
00493   if (pwc==0) {
00494     //get a string of the reweightVersion
00495     string sReweightVersion=this->AsString
00496       (static_cast<SKZPWeightCalculator::SKZPConfig_t>
00497        (nu.reweightVersion));
00498     MSG("NuZBeamReweight",Msg::kInfo)
00499       <<"To access database, using the string="<<sReweightVersion<<endl;
00500 
00501     TDirectory* tmpd = gDirectory;
00502     MSG("NuZBeamReweight",Msg::kDebug)
00503       <<"TDirectory before Reweight is:"<<endl;
00504     //tmpd->Print();
00505     
00506     //create the weight calculator
00507     //read constants from DB
00508     pwc=new SKZPWeightCalculator(sReweightVersion,true);
00509 
00510     //print the config
00511     pwc->PrintReweightConfig(cout);
00512 
00513     MSG("NuZBeamReweight",Msg::kDebug)
00514       <<"TDirectory after building Zbeam and Zflux is:"<<endl;
00515     //gDirectory->Print();
00516     MSG("NuZBeamReweight",Msg::kDebug)
00517       <<"After reseting gDirectory its location is:"<<endl;
00518     gDirectory=tmpd;
00519     //gDirectory->Print();
00520   }
00521   
00522   return *pwc;
00523 }
00524 
00525 //......................................................................
00526 
00527 void NuZBeamReweight::GetBeamWeightsOnly(NuEvent& nu) const
00528 {
00529   if (nu.simFlag!=SimFlag::kMC || 
00530       nu.reweightVersion==SKZPWeightCalculator::kUnknown) {
00531     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00532       <<"Not doing ExtractZBeamReweight, reweightVersion="
00533       <<nu.reweightVersion<<", simFlag="<<nu.simFlag<<endl;
00534     
00535     //set the weights to 1 that would otherwise be set below
00536     nu.trkEnWeight=1;//no reweighting
00537     nu.shwEnWeight=1;//no reweighting
00538     nu.beamWeight=1;//no reweighting
00539     nu.fluxErrHadProdAfterTune=-1;
00540     nu.fluxErrTotalErrorPreTune=-1;
00541     nu.fluxErrTotalErrorAfterTune=-1;
00542     nu.detectorWeightNMB=1;//no reweighting
00543     nu.detectorWeightNM=1;//no reweighting
00544     
00545     nu.trkEnWeightRunI=1;//no reweighting
00546     nu.shwEnWeightRunI=1;//no reweighting
00547     nu.beamWeightRunI=1;//no reweighting
00548     nu.fluxErrHadProdAfterTuneRunI=-1;
00549     nu.fluxErrTotalErrorPreTuneRunI=-1;
00550     nu.fluxErrTotalErrorAfterTuneRunI=-1;
00551     nu.detectorWeightNMBRunI=1;//no reweighting
00552     nu.detectorWeightNMRunI=1;//no reweighting
00553     
00554     nu.trkEnWeightRunII=1;//no reweighting
00555     nu.shwEnWeightRunII=1;//no reweighting
00556     nu.beamWeightRunII=1;//no reweighting
00557     nu.fluxErrHadProdAfterTuneRunII=-1;
00558     nu.fluxErrTotalErrorPreTuneRunII=-1;
00559     nu.fluxErrTotalErrorAfterTuneRunII=-1;
00560     nu.detectorWeightNMBRunII=1;//no reweighting
00561     nu.detectorWeightNMRunII=1;//no reweighting
00562 
00563     return;
00564   }
00565 
00566   MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00567     <<"Running GetZBeamReweight..."<<endl;
00568 
00569   //How many sigma for the beam errors?
00570   Float_t nsigma=1;
00571 
00572   //convert the std BeamType into Zarko's convention
00573   //was 2=le010z185i (Zbeam.h)
00574   const Int_t Ibeam=BeamType::ToZarko
00575     (static_cast<BeamType::BeamType_t>(nu.beamType));
00576     
00577   //get a reference
00578   SKZPWeightCalculator& wc=this->GetSKZPWeightCalculator(nu);
00579     
00580   //copied from SKZPWeightCalculator
00581   Int_t inu=nu.inu; 
00582   Double_t true_enu=nu.neuEnMC;//p4neu[3]
00583   Int_t parenttype = nu.tptype;
00584   Double_t parentpz = nu.tpz;
00585   Double_t parentpy = nu.tpy;
00586   Double_t parentpx = nu.tpx;
00587   Double_t pt = sqrt(parentpy*parentpy+parentpx*parentpx);
00588 
00590   //calc the weights for runI
00592   SKZPWeightCalculator::RunPeriod_t rp=SKZPWeightCalculator::kRunI;
00593     
00594   //get the skzp weights (beam only)
00595   Double_t beamweight1 = wc.GetBeamWeight(nu.detector,
00596                                           Ibeam,
00597                                           parenttype,
00598                                           pt,
00599                                           parentpz,
00600                                           true_enu,
00601                                           inu,
00602                                           rp);
00603 
00604   //get the flux errors
00605   Float_t fluxErrHadProdAfterTune1=wc.GetFluxError
00606     (nu.detector,Ibeam,inu,true_enu,
00607      SKZPWeightCalculator::kHadProdAfterTune,nsigma);
00608   Float_t fluxErrTotalErrorPreTune1=wc.GetFluxError
00609     (nu.detector,Ibeam,inu,true_enu,
00610      SKZPWeightCalculator::kTotalErrorPreTune,nsigma);
00611   Float_t fluxErrTotalErrorAfterTune1=wc.GetFluxError
00612     (nu.detector,Ibeam,inu,true_enu,
00613      SKZPWeightCalculator::kTotalErrorAfterTune,nsigma);
00614 
00616   //calc the weights for runII
00618   rp=SKZPWeightCalculator::kRunII;
00619     
00620   //get the skzp weights (beam only)
00621   Double_t beamweight2 = wc.GetBeamWeight(nu.detector,
00622                                           Ibeam,
00623                                           parenttype,
00624                                           pt,
00625                                           parentpz,
00626                                           true_enu,
00627                                           inu,
00628                                           rp);
00629 
00630   //get the flux errors
00631   Float_t fluxErrHadProdAfterTune2=wc.GetFluxError
00632     (nu.detector,Ibeam,inu,true_enu,
00633      SKZPWeightCalculator::kHadProdAfterTune,nsigma);
00634   Float_t fluxErrTotalErrorPreTune2=wc.GetFluxError
00635     (nu.detector,Ibeam,inu,true_enu,
00636      SKZPWeightCalculator::kTotalErrorPreTune,nsigma);
00637   Float_t fluxErrTotalErrorAfterTune2=wc.GetFluxError
00638     (nu.detector,Ibeam,inu,true_enu,
00639      SKZPWeightCalculator::kTotalErrorAfterTune,nsigma);
00640 
00642   //calc the weights for the weighted average
00644   static Bool_t firstTime=true;
00645   if (firstTime) {
00646     firstTime=false;
00647     //define and set the numbers of POTs in each run period
00648     //numbers from Andy Blakes minos-doc-3727
00649     Float_t potRunI=1.269e20;
00650     //runIIa=1.229
00651     //runIIb=0.721
00652     Float_t potRunII=1.950e20;
00653     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00654       <<"ExtractZBeamReweight: SetRunFrac w/ potRunI="<<potRunI
00655       <<", potRunII="<<potRunII<<endl;
00656     wc.SetRunFrac(Ibeam,SKZPWeightCalculator::kRunI,potRunI);
00657     wc.SetRunFrac(Ibeam,SKZPWeightCalculator::kRunII,potRunII);
00658   }
00659   
00660   //Beam weight will go in here
00661   Double_t beamweightAv = 0;
00662 
00663   //Make dummy variables to satisfy the detector weight part of the
00664   //code which we don't care about
00665   Int_t ccnc = 0;
00666   Double_t reco_emu = 0;
00667   Double_t reco_eshw = 0;
00668   Double_t new_reco_emu = 0;
00669   Double_t new_reco_eshw = 0;
00670   Double_t detweightAv = 0;
00671 
00672   //get the skzp weights (beam and detector)
00673   wc.GetRFWWeight(nu.detector,
00674                   Ibeam,
00675                   parenttype,
00676                   pt,
00677                   parentpz,
00678                   ccnc,
00679                   true_enu,
00680                   inu,
00681                   reco_emu,
00682                   reco_eshw,
00683                   new_reco_emu,
00684                   new_reco_eshw,
00685                   beamweightAv,
00686                   detweightAv);
00687   
00688   //get the flux errors
00689   Float_t fluxErrHadProdAfterTune=wc.GetFluxError
00690     (nu.detector,Ibeam,inu,true_enu,
00691      SKZPWeightCalculator::kHadProdAfterTune,nsigma);
00692   Float_t fluxErrTotalErrorPreTune=wc.GetFluxError
00693     (nu.detector,Ibeam,inu,true_enu,
00694      SKZPWeightCalculator::kTotalErrorPreTune,nsigma);
00695   Float_t fluxErrTotalErrorAfterTune=wc.GetFluxError
00696     (nu.detector,Ibeam,inu,true_enu,
00697      SKZPWeightCalculator::kTotalErrorAfterTune,nsigma);
00698 
00699   MAXMSG("NuZBeamReweight",Msg::kDebug,300)
00700     <<"Wgts RunI:  b="<<beamweight1 << endl
00701     <<"Wgts RunII: b="<<beamweight2 << endl
00702     <<"Wgts Av:    b="<<beamweightAv << endl;
00703   
00704   MAXMSG("NuZBeamReweight",Msg::kDebug,200)
00705     <<"energyMC="<<nu.neuEnMC
00706     <<", inu="<<nu.inu
00707     <<endl
00708     <<"fluxErrs:"<<" HadProdAfterTune1   ="
00709     <<fluxErrHadProdAfterTune1<<endl
00710     <<"         "<<" TotalErrorPreTune1  ="
00711     <<fluxErrTotalErrorPreTune1<<endl
00712     <<"         "<<" TotalErrorAfterTune1="
00713     <<fluxErrTotalErrorAfterTune1<<endl
00714     <<"fluxErrs:"<<" HadProdAfterTune2   ="
00715     <<fluxErrHadProdAfterTune2<<endl
00716     <<"         "<<" TotalErrorPreTune2  ="
00717     <<fluxErrTotalErrorPreTune2<<endl
00718     <<"         "<<" TotalErrorAfterTune2="
00719     <<fluxErrTotalErrorAfterTune2<<endl
00720     <<"fluxErrs:"<<" HadProdAfterTune   ="
00721     <<fluxErrHadProdAfterTune<<endl
00722     <<"         "<<" TotalErrorPreTune  ="
00723     <<fluxErrTotalErrorPreTune<<endl
00724     <<"         "<<" TotalErrorAfterTune="
00725     <<fluxErrTotalErrorAfterTune<<endl;
00726  
00728   //store the weights
00730   nu.beamWeight=beamweightAv;
00731   nu.fluxErrHadProdAfterTune=fluxErrHadProdAfterTune;
00732   nu.fluxErrTotalErrorPreTune=fluxErrTotalErrorPreTune;
00733   nu.fluxErrTotalErrorAfterTune=fluxErrTotalErrorAfterTune;
00734     
00735   nu.beamWeightRunI=beamweight1;
00736   nu.fluxErrHadProdAfterTuneRunI=fluxErrHadProdAfterTune1;
00737   nu.fluxErrTotalErrorPreTuneRunI=fluxErrTotalErrorPreTune1;
00738   nu.fluxErrTotalErrorAfterTuneRunI=fluxErrTotalErrorAfterTune1;
00739     
00740   nu.beamWeightRunII=beamweight2;
00741   nu.fluxErrHadProdAfterTuneRunII=fluxErrHadProdAfterTune2;
00742   nu.fluxErrTotalErrorPreTuneRunII=fluxErrTotalErrorPreTune2;
00743   nu.fluxErrTotalErrorAfterTuneRunII=fluxErrTotalErrorAfterTune2;
00744 }
00745 
00746 //......................................................................
00747 
00748 void NuZBeamReweight::CalcGeneratorReweight(const NtpStRecord& ntp,
00749                                             const NtpSREvent& evt,
00750                                             NuEvent& nu) const
00751 {
00755 
00756   //return if not MC
00757   if (nu.simFlag!=SimFlag::kMC) {
00758     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00759       <<"Not doing CalcGeneratorReweight"
00760       <<"for simFlag="<<nu.simFlag<<endl;
00761     
00762     //set the weights to 1 that would otherwise be set below
00763     nu.generatorWeight=1;//no reweighting
00764     return;
00765   }
00766  
00767   if (!nu.useGeneratorReweight) {
00768     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00769       <<"CalcGeneratorReweight: not calculating Generator reweight"
00770       <<" since nu.useGeneratorReweight="<<nu.useGeneratorReweight
00771       <<endl;
00772     return;
00773   }
00774 
00775   if (ReleaseType::IsDaikon(nu.releaseType)) {
00776     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00777       <<"Not calculating Generator weight for Daikon, not required"
00778       <<endl;
00779     return;
00780   }
00781   
00782   MAXMSG("NuZBeamReweight",Msg::kDebug,100)
00783     <<"CalcGeneratorReweight:"<<endl
00784     <<"  generator:config_name="<<nu.sGeneratorConfigName
00785     <<"  generator:config_no="<<nu.generatorConfigNo
00786     <<endl;
00787   
00788   //declare the registry and create it only once
00789   static Registry* rwtconfig=this->CreateNeugenRegistry(nu);
00790 
00791   //make a WeightCalculator and add it to the MCReweight singleton
00792   this->AddNeugenWeightCalculator();
00793 
00794   //get the MC info from the ntp
00795   TClonesArray& thevtTca=*ntp.thevt;
00796   const NtpTHEvent* the=
00797     dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
00798   Int_t thidx=the->neumc;
00799 
00800   //fill MCEventInfo to do MODBYRS3;
00801   MCEventInfo ei;
00802   ei.UseStoredXSec(true);
00803   NtpStRecord* st=const_cast<NtpStRecord*>(&ntp);//nasty!
00804   ReweightHelpers::MCEventInfoFilla(&ei,st,thidx);
00805 
00806   //a second one for a sanity check on the values in the NuEvent
00807   MCEventInfo ei2;
00808   ei2.UseStoredXSec(true);
00809   this->MCEventInfoFilla(&ei2,nu);
00810   this->CompareMCEventInfo(ei,ei2);
00811   
00812   //Compute the MODBYRS3 weight
00813   Double_t generatorweight=1.;
00814   NuParent* np=0;
00815   if (nu.useGeneratorReweight){
00816     MCReweight* mcr=&MCReweight::Instance();
00817     generatorweight=mcr->ComputeWeight(&ei,np,rwtconfig);
00818   }
00819   
00820   //set the output of the function
00821   nu.generatorWeight=generatorweight;
00822 }
00823 
00824 //......................................................................
00825 
00826 void NuZBeamReweight::CalcGeneratorReweight(NuEvent& nu) const
00827 {
00828   if (!nu.useGeneratorReweight) {
00829     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00830       <<"CalcGeneratorReweight(nu): not calculating Generator reweight"
00831       <<" since nu.useGeneratorReweight="<<nu.useGeneratorReweight
00832       <<endl;
00833     return;
00834   }
00835 
00836   if (ReleaseType::IsDaikon(nu.releaseType)) {
00837     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00838       <<"Not doing CalcGeneratorReweight(nu) for Daikon, not required"
00839       <<endl;
00840     return;
00841   }
00842 
00843   MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00844     <<"CalcGeneratorReweight(nu):"<<endl
00845     <<"  generator:config_name="<<nu.sGeneratorConfigName
00846     <<"  generator:config_no="<<nu.generatorConfigNo
00847     <<endl;
00848   
00849   //declare the registry and create it only once
00850   static Registry* rwtconfig=this->CreateNeugenRegistry(nu);
00851 
00852   //make a WeightCalculator and add it to the MCReweight singleton
00853   this->AddNeugenWeightCalculator();
00854   
00855   //fill MCEventInfo
00856   MCEventInfo ei;
00857   ei.UseStoredXSec(true);
00858   //ReweightHelpers::MCEventInfoFilla(&ei,st,thidx);
00859   this->MCEventInfoFilla(&ei,nu);
00860   
00861   //Compute the MODBYRS3 weight
00862   Double_t generatorweight=1.;
00863   NuParent* np=0;
00864   if (nu.useGeneratorReweight){
00865     MCReweight* mcr=&MCReweight::Instance();
00866     generatorweight=mcr->ComputeWeight(&ei,np,rwtconfig);
00867   }
00868   
00869   //set the output of the function
00870   nu.generatorWeight=generatorweight;
00871 }
00872 
00873 //......................................................................
00874 
00875 void NuZBeamReweight::AddNeugenWeightCalculator() const
00876 {
00878   Bool_t firstTime=true;
00879   if (firstTime) {
00880     firstTime=false;
00881     MAXMSG("NuZBeamReweight",Msg::kInfo,1)
00882       <<"Adding NeugenWeightCalculator to MCReweight singleton"<<endl;
00883     static NeugenWeightCalculator *n=new NeugenWeightCalculator();
00884     MCReweight* mcr=&MCReweight::Instance();    
00885     mcr->AddWeightCalculator(n);
00886   }
00887 }
00888 
00889 //......................................................................
00890 
00891 Registry* NuZBeamReweight::CreateNeugenRegistry(const NuEvent& nu) const
00892 {
00893   static Registry* rwtconfig=0;
00894   
00895   //create it the first time
00896   if (rwtconfig==0){
00897     
00898     //fill the registry to tell MCReweight to do MODBYRS3
00899     rwtconfig=new Registry();
00900     rwtconfig->UnLockValues();
00901     rwtconfig->UnLockKeys();
00902     rwtconfig->Clear();
00903     rwtconfig->Set("neugen:config_name",nu.sGeneratorConfigName.c_str());
00904     rwtconfig->Set("neugen:config_no",nu.generatorConfigNo);
00905     rwtconfig->LockValues();
00906     rwtconfig->LockKeys();
00907   }
00908   return rwtconfig;
00909 }
00910 
00911 //......................................................................
00912 
00913 void NuZBeamReweight::CompareMCEventInfo(const MCEventInfo& ei,
00914                                          const MCEventInfo& ei2) const
00915 {
00916   if (ei.initial_state!=ei2.initial_state ||
00917       ei.nucleus!=ei2.nucleus ||
00918       ei.had_fs!=ei2.had_fs){
00919     MAXMSG("NuZBeamReweight",Msg::kError,300)
00920       <<"Ahhh, CompareMCEventInfo difference found:"<<endl
00921       <<"1st: initial_state="<<ei.initial_state
00922       <<", nucleus="<<ei.nucleus
00923       <<", had_fs="<<ei.had_fs
00924       <<endl
00925       <<"2nd: initial_state="<<ei2.initial_state
00926       <<", nucleus="<<ei2.nucleus
00927       <<", had_fs="<<ei2.had_fs
00928       <<endl
00929       <<"NOTE: This discrepancy has to be understood if we"
00930       <<" are to trust the generator reweighting"
00931       <<endl;
00932   }
00933 }
00934 
00935 //......................................................................
00936 
00937 void NuZBeamReweight::MCEventInfoFilla(MCEventInfo* ei,
00938                                        const NuEvent& nu) const
00939 {
00942 
00943   ei->nuE=nu.neuEnMC;
00944   ei->nuPx=nu.neuPxMC;
00945   ei->nuPy=nu.neuPyMC;
00946   ei->nuPz=nu.neuPzMC;
00947   
00948   ei->tarE=nu.tgtEnMC;
00949   ei->tarPx=nu.tgtPxMC;
00950   ei->tarPy=nu.tgtPyMC;
00951   ei->tarPz=nu.tgtPzMC;
00952   
00953   ei->y=nu.yMC;
00954   ei->x=nu.xMC;
00955   ei->q2=nu.q2MC;
00956   ei->w2=nu.w2MC;
00957   
00958   ei->iaction=nu.iaction;
00959   ei->inu=nu.inu;
00960   ei->iresonance=nu.iresonance;
00961   ei->initial_state=nu.initialStateMC;
00962   
00963   ei->nucleus=nu.nucleusMC;
00964   
00965   ei->had_fs=nu.hadronicFinalStateMC;
00966 
00967 
00968   MAXMSG("NuZBeamReweight",Msg::kDebug,100)
00969     <<endl<<endl
00970     <<", ei->nuE="<<ei->nuE
00971     <<", ei->nuPx="<<ei->nuPx
00972     <<", ei->nuPy="<<ei->nuPy
00973     <<", ei->nuPz="<<ei->nuPz
00974     <<endl
00975     <<", ei->tarE="<<ei->tarE
00976     <<", ei->tarPx="<<ei->tarPx
00977     <<", ei->tarPy="<<ei->tarPy
00978     <<", ei->tarPz="<<ei->tarPz
00979     <<endl
00980     <<", ei->y="<<ei->y<<endl
00981     <<", ei->x="<<ei->x<<endl
00982     <<", ei->q2="<<ei->q2<<endl
00983     <<", ei->w2="<<ei->w2<<endl
00984     <<endl
00985     <<", ei->iaction="<<ei->iaction<<endl
00986     <<", ei->inu="<<ei->inu<<endl
00987     <<", ei->iresonance="<<ei->iresonance<<endl
00988     <<", ei->initial_state="<<ei->initial_state<<endl
00989     <<", ei->nucleus="<<ei->nucleus<<endl
00990     <<", ei->had_fs="<<ei->had_fs<<endl
00991     <<endl;
00992 }
00993 
00994 //......................................................................
00995 
00996 const Char_t* NuZBeamReweight::
00997 AsString(SKZPWeightCalculator::SKZPConfig_t c) const
00998 {
00999   // return char string in cannonical offline form
01000   switch (c) {
01001     case SKZPWeightCalculator::kPRL:     return "PRL";
01002     case SKZPWeightCalculator::kBoston:  return "Boston";
01003     case SKZPWeightCalculator::kPiMinus: return "PiMinus_CedarDaikon";
01004     case SKZPWeightCalculator::kDetXs:   return "DetXs";
01005     case SKZPWeightCalculator::kUnknown: return "Unknown"; break;
01006     default:         return "!?Bad Value Passed?!"; break;
01007   }
01008 }
01009 
01010 //......................................................................