NueReweight Class Reference

#include <NueReweight.h>

Inheritance diagram for NueReweight:
JobCModule

List of all members.

Public Member Functions

 NueReweight ()
 ~NueReweight ()
void EndJob ()
int ReadRandom ()
virtual void Reset ()
JobCResult Reco (MomNavigator *mom)
const RegistryDefaultConfig () const
void Config (const Registry &r)

Public Attributes

NueRWrw

Private Attributes

MCReweight mcr
int counter
int rand
bool firstevent
Registry rwtconfig

Detailed Description

Definition at line 19 of file NueReweight.h.


Constructor & Destructor Documentation

NueReweight::NueReweight (  ) 

Definition at line 31 of file NueReweight.cxx.

References Msg::kDebug, MSG, and n.

00031                         :
00032    mcr(MCReweight::Instance()),
00033    counter(0),
00034    rand(0),
00035    firstevent(true)
00036 {
00037    MSG("NueReweight",Msg::kDebug)<<"Making a NueReweight"<<endl;
00038    rw = new NueRW();
00039 
00040    // make a WeightCalculator and add it to the MCReweight singleton
00041    NeugenWeightCalculator *n=new NeugenWeightCalculator();
00042    mcr.AddWeightCalculator(n);
00043 }

NueReweight::~NueReweight (  ) 

Definition at line 47 of file NueReweight.cxx.

References Msg::kDebug, MSG, and rw.

00048 {
00049    MSG("NueReweight",Msg::kDebug)<<"Killing a NueReweight"<<endl;
00050    if(rw){
00051       delete rw;
00052       rw=0;
00053    }
00054 }


Member Function Documentation

void NueReweight::Config ( const Registry r  )  [virtual]

Return the actual configuration. If your module directly pulls its configuration from the fConfig Registry, you don't need to override this. Override if you have local config variables.

Reimplemented from JobCModule.

Definition at line 247 of file NueReweight.cxx.

References Registry::Get(), and rand.

00248 {
00249   int    tmpi;
00250 
00251   if (r.Get("RandNumber",tmpi)) { rand = tmpi; }
00252 }

const Registry & NueReweight::DefaultConfig ( void   )  const [virtual]

Get the default configuration registry. This should normally be overridden. One useful idiom is to implement it like:

const Registry& MyModule::DefaultConfig() const { static Registry cfg; // never is destroyed if (cfg.Size()) return cfg; // already filled it // set defaults: cfg.Set("TheAnswer",42); cfg.Set("Units","unknown"); return cfg; }

Reimplemented from JobCModule.

Definition at line 228 of file NueReweight.cxx.

References JobCModule::GetName(), Registry::LockValues(), Registry::Set(), and Registry::UnLockValues().

00229 {
00230   static Registry r; // Default configuration for module
00231 
00232   // Set name of config
00233   std::string name = this->GetName();
00234   name += ".config.default";
00235   r.SetName(name.c_str());
00236 
00237   // Set values in configuration
00238   r.UnLockValues();
00239   r.Set("RandNumber", 0);
00240   r.LockValues();
00241 
00242   return r;
00243 }

void NueReweight::EndJob (  )  [virtual]

Implement for notification of end of job

Reimplemented from JobCModule.

Definition at line 58 of file NueReweight.cxx.

References Msg::kDebug, and MSG.

00059 {
00060    MSG("NueReweight",Msg::kDebug)<<"End Job"<<endl;
00061 }

int NueReweight::ReadRandom (  ) 

Definition at line 257 of file NueReweight.cxx.

References NueRW::coh_ma, NueRW::coh_r0, NueRW::coh_rei, counter, NueRW::dm2, fname, Msg::kError, NueRW::kno_a1, NueRW::kno_a2, NueRW::kno_a3, NueRW::kno_a4, NueRW::kno_b, NueRW::kno_r112, NueRW::kno_r113, NueRW::kno_r122, NueRW::kno_r123, NueRW::kno_r132, NueRW::kno_r133, NueRW::kno_r142, NueRW::kno_r143, MSG, NueRW::qel_eta, NueRW::qel_fa0, NueRW::qel_ma, NueRW::randrow, NueRW::res_ma, NueRW::res_omega, NueRW::res_z, rw, NueRW::ss2th, and NueRW::UE32.

Referenced by Reco().

00258 {
00259    const Int_t NCOLS=25;
00260 
00261    string fname(getenv("SRT_PRIVATE_CONTEXT"));
00262    //   fname+=("/NueAna/data/randnumbers.txt");
00263    //   fname+=("/NueAna/data/uniformnumbers.txt");
00264    fname+=("/NueAna/data/biguniformnumbers.txt");
00265    ifstream in(fname.c_str());
00266    if(!in){
00267       MSG("NueReweight",Msg::kError)<<"Could not open random number file in priviate area"<<endl
00268                                     <<"Tried: "<<fname<<endl;
00269       fname=getenv("SRT_PUBLIC_CONTEXT");
00270       //      fname+=("/NueReweight/data/randnumbers.txt");
00271       fname+=("/NueAna/data/biguniformnumbers.txt");
00272       MSG("NueReweight",Msg::kError)<<"Will now try: "<<fname<<endl;
00273       in.open(fname.c_str());
00274       if(!in){
00275          MSG("NueReweight",Msg::kError)<<"Could not open a random number file"<<endl;
00276          return -1;
00277       }
00278    }
00279 
00280    float rands[NCOLS];
00281    int counter=0;
00282    in>>rands[0];
00283    while(!in.eof()){
00284       for(int i=1;i<NCOLS;i++){
00285          in>>rands[i];
00286       }
00287       if(counter==rw->randrow){
00288          break;
00289       }
00290       counter++;
00291       in>>rands[0];
00292    }
00293 
00294    if(counter<rw->randrow){
00295       MSG("NueReweight",Msg::kError)<<"Couldn't get row "<<rw->randrow<<" from "<<fname<<endl;
00296       return -1;
00297    }
00298    
00299    rw->qel_ma=fabs(rands[0]);
00300    rw->res_ma=fabs(rands[1]);
00301    rw->coh_ma=fabs(rands[2]);
00302    rw->qel_fa0=-1.*fabs(rands[3]);
00303    rw->qel_eta=fabs(rands[4]);
00304    rw->res_omega=fabs(rands[5]);
00305    rw->res_z=fabs(rands[6]);
00306    rw->coh_r0=fabs(rands[7]);
00307    rw->coh_rei=fabs(rands[8]);
00308    rw->kno_a1=fabs(rands[9]);
00309    rw->kno_a2=fabs(rands[10]);
00310    rw->kno_a3=fabs(rands[11]);
00311    rw->kno_a4=fabs(rands[12]);
00312    rw->kno_b=fabs(rands[13]);
00313 
00314    //use same scale for all kno_r* variables
00315    rw->kno_r112=fabs(rands[14]);
00316    rw->kno_r122=fabs(rands[14]);
00317    rw->kno_r132=fabs(rands[14]);
00318    rw->kno_r142=fabs(rands[14]);   
00319    rw->kno_r113=fabs(rands[14]);
00320    rw->kno_r123=fabs(rands[14]);
00321    rw->kno_r133=fabs(rands[14]);
00322    rw->kno_r143=fabs(rands[14]);
00323    rw->dm2=fabs(rands[22]);
00324    if(rands[23]>1){
00325      rw->ss2th=1.-(rands[23]-1.);
00326    }
00327    else{
00328      rw->ss2th=fabs(rands[23]);
00329    }
00330    rw->UE32=fabs(rands[24]);
00331    return 0;
00332 }

JobCResult NueReweight::Reco ( MomNavigator mom  )  [virtual]

Implement this for read-write access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 65 of file NueReweight.cxx.

References MCReweight::ComputeWeight(), counter, NueRW::dm2, ANtpEventInfo::energyGeV, ReweightHelpers::EventRegistryFilla(), NueRW::fDet, NueRW::fFileType, NueRW::FindEBin(), firstevent, MomNavigator::FragmentIter(), NueRW::fRun, NueRW::fSubRun, VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), NueHeader::GetRun(), NueHeader::GetSubRun(), RecHeader::GetVldContext(), ANtpTruthInfo::interactionType, NueRW::kBEAM, Msg::kDebug, Detector::kFar, NueRW::kno_r112, NueRW::kno_r113, NueRW::kno_r122, NueRW::kno_r123, NueRW::kno_r132, NueRW::kno_r133, NueRW::kno_r142, NueRW::kno_r143, NueRW::kNUE, JobCResult::kPassed, NueRW::kTAU, NueRW::kUnknown, Registry::LockKeys(), Registry::LockValues(), mcr, NueRecord::mctrue, MSG, NueRW::nacc, NueRW::nbg, NueRW::nbgE, NueRW::nevents, NueRW::nnc, NueRW::nncE, NueRW::nnueb, NueRW::nnuebE, NueRW::nnumu, NueRW::nnumuE, NueRW::nnutau, NueRW::nnutauE, ANtpTruthInfoBeam::nonOscNuFlavor, NueRW::nsig, NueRW::nsigE, NueRW::nsnarls, ANtpTruthInfo::nuFlavor, Oscillate(), NueRW::qel_ma, rand, NueRW::randrow, ReadRandom(), NueRW::res_ma, NueRW::Reset(), rw, rwtconfig, Registry::Set(), NueRecord::srevent, NueRW::ss2th, NueRW::UE32, Registry::UnLockKeys(), and Registry::UnLockValues().

00066 {
00067   //  if(counter%1000==0){
00068   //      MSG("NueReweight",Msg::kInfo)<<"On entry "<<counter<<endl;
00069   //  }
00070   rw->randrow=rand;
00071    counter++;
00072    TObject *obj=0;
00073    TIter objiter = mom->FragmentIter();
00074    int objcount=0;
00075    while((obj=objiter.Next())){
00076       const char *cn=obj->ClassName();
00077       //      MSG("NueReweight",Msg::kDebug)<<"Found a "<<cn
00078       //                            <<". Objcount "<<objcount<<endl;
00079       if(strcmp(cn,"NueRecord")!=0){
00080         continue;
00081       }
00082       objcount++;
00083       NueRecord *nr = dynamic_cast<NueRecord *>(obj);
00084       if(firstevent){
00085         MSG("NueReweight",Msg::kDebug)<<"this is the first event"<<endl;
00086         rw->Reset();
00087         rw->fRun = nr->GetHeader().GetRun();
00088         rw->fSubRun = nr->GetHeader().GetSubRun();
00089         rw->fDet=nr->GetHeader().GetVldContext().GetDetector();
00090         int iflv = (int)((rw->fRun%100000-rw->fRun%1000)/10000);
00091         MSG("NueReweight",Msg::kDebug)<<"run is "<<rw->fRun       
00092                                       <<" run%100000 is "<<rw->fRun%100000
00093                                       <<" run%1000 is "<<rw->fRun%1000<<endl;
00094         
00095         MSG("NueReweight",Msg::kDebug)<<"iflv is "<<iflv<<endl;
00096         if(iflv==0){
00097           rw->fFileType = NueRW::kBEAM;
00098         }
00099         else if(iflv==1){
00100           rw->fFileType = NueRW::kNUE;
00101         }
00102         else if(iflv==3){
00103           rw->fFileType = NueRW::kTAU;
00104         }
00105         else{
00106           rw->fFileType = NueRW::kUnknown;
00107         }
00108         rw->randrow=rand;
00109         ReadRandom();      
00110       // to calculate a reweight factor, two registry objects are used
00111       // the first contains the list of parameters in the WeightCalculator to
00112       // vary 
00113       // the list of parameters that can be varied can be seen by looking in the
00114       // appropriate WeightCalculator::Config() methods
00115 
00116         rwtconfig.UnLockValues();
00117         rwtconfig.UnLockKeys();
00118         rwtconfig.Set("neugen:use_scale_factors",1);
00119         rwtconfig.Set("neugen:ma_qe",rw->qel_ma);
00120         rwtconfig.Set("neugen:ma_res",rw->res_ma);
00121         //      rwtconfig.Set("neugen:ma_coh",rw->coh_ma);
00122         //      rwtconfig.Set("neugen:qel_fa0",rw->qel_fa0);
00123         //      rwtconfig.Set("neugen:qek_eta",rw->qel_eta);
00124         //      rwtconfig.Set("neugen:res_qe",rw->res_omega);
00125         //      rwtconfig.Set("neugen:res_z",rw->res_z);
00126         //      rwtconfig.Set("neugen:coh_r0",rw->coh_r0);
00127         //      rwtconfig.Set("neugen:coh_rei",rw->coh_rei);
00128         //      rwtconfig.Set("neugen:kno_a1",rw->kno_a1);
00129         //      rwtconfig.Set("neugen:kno_a2",rw->kno_a2);
00130         //      rwtconfig.Set("neugen:kno_a3",rw->kno_a3);
00131         //      rwtconfig.Set("neugen:kno_a4",rw->kno_a4);
00132         //      rwtconfig.Set("neugen:kno_b",rw->kno_b);
00133         rwtconfig.Set("neugen:kno_r112",rw->kno_r112);
00134         rwtconfig.Set("neugen:kno_r122",rw->kno_r122);
00135         rwtconfig.Set("neugen:kno_r132",rw->kno_r132);
00136         rwtconfig.Set("neugen:kno_r142",rw->kno_r142);
00137         rwtconfig.Set("neugen:kno_r113",rw->kno_r113);
00138         rwtconfig.Set("neugen:kno_r123",rw->kno_r123);
00139         rwtconfig.Set("neugen:kno_r133",rw->kno_r133);
00140         rwtconfig.Set("neugen:kno_r143",rw->kno_r143);
00141         rwtconfig.LockValues();
00142         rwtconfig.LockKeys();
00143 
00144         firstevent=false;
00145         MSG("NueReweight",Msg::kDebug)<<"Set rw->randrow to "<<rw->randrow<<endl;
00146       }
00147       if(objcount==1){
00148         rw->nsnarls++;
00149       }
00150       rw->nevents++;   
00151       rw->nacc++;
00152       //do reweighting--
00153 
00154       //make event regsitry
00155       Registry evreg;
00156       ReweightHelpers::EventRegistryFilla(&(nr->mctrue),evreg);
00157 
00158       // now all that's left is to compute the weight as follows:
00159       float reweight = mcr.ComputeWeight(&evreg,&rwtconfig);
00160 
00161       MSG("NueReweight",Msg::kDebug)<<"Reweight is "<<reweight<<endl;
00162        
00163       MSG("NueReweight",Msg::kDebug)<<"*********************"<<endl;      
00164       //      evreg.Print();
00165       MSG("NueReweight",Msg::kDebug)<<"*********************"<<endl;      
00166       //      rwtconfig.Print();
00167 
00168 
00169       //oscillate
00170       float oscprob = 1.;
00171       if(nr->GetHeader().GetVldContext().GetDetector()==Detector::kFar){
00172         float theta23 = TMath::ASin(sqrt(rw->ss2th))/2.;
00173         oscprob = NueConvention::Oscillate(&(nr->mctrue),735.,
00174                                           rw->dm2,theta23,rw->UE32);
00175         MSG("NueReweight",Msg::kDebug)<<"dm2 is "<<rw->dm2<<" theta23 is "<<theta23
00176                                       <<" ue32 is "<<rw->UE32<<endl;
00177       }
00178       MSG("NueReweight",Msg::kDebug)<<"oscprob is "<<oscprob<<endl;
00179       int ebin=rw->FindEBin(nr->srevent.energyGeV);
00180       float scale=reweight*oscprob;
00181 
00182       int iaction = nr->mctrue.interactionType;
00183       int inu = nr->mctrue.nuFlavor;
00184       int inunoosc = nr->mctrue.nonOscNuFlavor;
00185       
00186       if(iaction==0){
00187          //neutral current;
00188         rw->nnc+=scale;
00189         rw->nbg+=scale;
00190         rw->nncE[ebin]+=scale;
00191         rw->nbgE[ebin]+=scale;
00192       }
00193       else{
00194         if(abs(inu)==12){
00195           if(abs(inunoosc)==12){
00196              rw->nnueb+=scale;
00197              rw->nbg+=scale;
00198              rw->nnuebE[ebin]+=scale;
00199              rw->nbgE[ebin]+=scale;
00200           }
00201           else if(abs(inunoosc)==14){
00202             rw->nsig+=scale;
00203             rw->nsigE[ebin]+=scale;
00204           }
00205         }
00206         else if(abs(inu)==14){
00207           rw->nnumu+=scale;
00208           rw->nbg+=scale;
00209           rw->nnumuE[ebin]+=scale;
00210           rw->nbgE[ebin]+=scale;
00211         }
00212         else if(abs(inu)==16){
00213           rw->nnutau+=scale;
00214           rw->nbg+=scale;
00215           rw->nnutauE[ebin]+=scale;
00216           rw->nbgE[ebin]+=scale;
00217         }
00218       }
00219    }          
00220 //   MSG("NueReweight",Msg::kDebug)<<"InReco"<<endl;
00221    
00222 
00223    return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00224 }

void NueReweight::Reset (  )  [virtual]

Implement to reset oneself

Reimplemented from JobCModule.

Definition at line 336 of file NueReweight.cxx.

References Registry::Clear(), counter, firstevent, Registry::LockKeys(), Registry::LockValues(), rand, NueRW::Reset(), rw, rwtconfig, Registry::UnLockKeys(), and Registry::UnLockValues().

00337 {
00338    counter=0;
00339    rand=0;
00340    firstevent=true;
00341    rw->Reset();
00342    rwtconfig.UnLockKeys();
00343    rwtconfig.UnLockValues();
00344    rwtconfig.Clear();
00345    rwtconfig.LockKeys();
00346    rwtconfig.LockValues();
00347    return;
00348 }


Member Data Documentation

int NueReweight::counter [private]

Definition at line 44 of file NueReweight.h.

Referenced by ReadRandom(), Reco(), and Reset().

bool NueReweight::firstevent [private]

Definition at line 46 of file NueReweight.h.

Referenced by Reco(), and Reset().

Definition at line 43 of file NueReweight.h.

Referenced by Reco().

int NueReweight::rand [private]

Definition at line 45 of file NueReweight.h.

Referenced by Config(), Reco(), and Reset().

Definition at line 40 of file NueReweight.h.

Referenced by ReadRandom(), Reco(), Reset(), and ~NueReweight().

Definition at line 47 of file NueReweight.h.

Referenced by Reco(), and Reset().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1