loon: NueFluxWeightsAna.cxx Source File

00001 #include "StandardNtuple/NtpStRecord.h"
00002 #include "NueAna/NueAnaTools/SntpHelpers.h"
00003 #include "MCNtuple/NtpMCTruth.h"
00004 #include "MessageService/MsgService.h"
00005 #include "MCReweight/Zbeam.h"
00006 #include "MCReweight/Zfluk.h"
00007 #include "MCReweight/Kfluk.h"
00008 #include "NueAna/NueFluxWeightsAna.h"
00009 
00010 CVSID("$Id: NueFluxWeightsAna.cxx,v 1.14 2009/06/24 22:43:52 vahle Exp $");
00011 
00012 NueFluxWeightsAna::NueFluxWeightsAna(NueFluxWeights &nuefw):
00013   fNueFluxWeight(nuefw),
00014 //  beam(2),
00015   det(1)
00016 {
00017   //   cfg = "PiMinus_CedarDaikon";
00018    cfg = "DetXs";
00019 }
00020 
00021 NueFluxWeightsAna::~NueFluxWeightsAna()
00022 {
00023   MSG("NueFluxWeightsAns",Msg::kDebug)<<"in NueFluxWeightsAna destructor"<<endl;
00024 }
00025 
00026 void NueFluxWeightsAna::Analyze(int evtn, RecRecordImp<RecCandHeader> *srobj)
00027 {
00028   NtpStRecord *st=dynamic_cast<NtpStRecord *>(srobj);
00029   if(st==0){
00030     MSG("NueFluxWeightsAna",Msg::kError)<<"Trying to do flux reweighting on an event"
00031                                         <<" that comes from a non NtpStRecord"<<endl
00032                                         <<"That's not a good idea"<<endl;
00033     return;
00034   }
00035 
00036   int thn = SntpHelpers::GetEvent2MCIndex(evtn,st);
00037   NtpMCTruth *mcrec = SntpHelpers::GetMCTruth(thn,st);
00038   if(mcrec==0){
00039     return;
00040   }
00041   MSG("NueFluxWeightsAns",Msg::kDebug)<<"in NueFluxWeightsAna::Analyze"<<endl;
00042 
00043   if(fi==0){
00044     MSG("NueFluxWeightsAna",Msg::kWarning)<<"No FluxInfo object set, "
00045                                           <<"using the one from NtpMCTruth"
00046                                           <<" Once the flux files are fixed, "
00047                                           <<"and we no longer need the MuPi trees, "
00048                                           <<"we can comment out this comment"<<endl;
00049     fi=&mcrec->flux;
00050   }
00051 
00052   MSG("NueFluxWeightsAna",Msg::kDebug)<<"starting flux reweight"<<endl;
00053 
00054   double pt = sqrt(fi->tpx*fi->tpx+fi->tpy*fi->tpy);
00055   double pz = 1.*fi->tpz;
00056   int tptype = fi->tptype;
00057   int inu = mcrec->inu; 
00058   int cc_nc = mcrec->iaction;
00059 
00060   float true_enu = mcrec->p4neu[3];
00061 
00062   zbeam = BeamType::ToZarko(beam);
00063   MSG("NueFluxWeightsAna",Msg::kDebug)<<"true_enu "<<true_enu<<" beam (Z) "<<zbeam<<" det "<<det<<endl;
00064   double kflukw = 1.;//kfluk->GetWeight(fi);
00065   fNueFluxWeight.kflukweight = kflukw;
00066 
00067   double newTrackE, newShwE;
00068   double bweight =  skzpCalc->GetBeamWeight(det,zbeam,tptype,pt,pz,true_enu,inu);
00069   double dweight =  skzpCalc->GetDetWeight(cc_nc,true_enu,inu,0,0,newTrackE,newShwE);
00070   
00071   fNueFluxWeight.totbeamweight = bweight;
00072   fNueFluxWeight.totskzpweight = dweight*bweight;
00073   fNueFluxWeight.detectorWeight = dweight;
00074 
00075   fNueFluxWeight.skzpTrkEnergy = 0.0;  //newTrackE;
00076   fNueFluxWeight.skzpShwEnergy = 0.0; //newShwE;
00077   fNueFluxWeight.skzpConfig = cfg;
00078 
00079   for(int rpit=SKZPWeightCalculator::kRunI;rpit<SKZPWeightCalculator::kEndOfList;rpit++){
00080     SKZPWeightCalculator::RunPeriod_t rp = skzpCalc->RunPeriodFromInt(rpit);
00081     double rpw = skzpCalc->GetBeamWeight(det,zbeam,tptype,pt,pz,true_enu,inu,rp);
00082     fNueFluxWeight.RPtotbeamweight.push_back(rpw);
00083   }
00084 
00085 
00086   MSG("NueFluxWeightsAna",Msg::kDebug)<<"alldone with flux reweight "<<fNueFluxWeight.totbeamweight<<endl;
00087 }
00088