loon: NueReadTJPID.cxx Source File

00001 
00002 // $Id: NueReadTJPID.cxx,v 1.1 2006/09/18 21:41:22 boehm Exp $
00003 //
00004 // FILL_IN: [Document your code!!]
00005 //
00006 // vahle@hep.ucl.ac.uk
00008 #include "TH1F.h"
00009 #include "TFile.h"
00010 #include "NueAna/NueReadTJPID.h"
00011 #include "NueAna/NuePID.h"
00012 #include "MessageService/MsgService.h"
00013 #include "MinosObjectMap/MomNavigator.h"
00014 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00015 #include "HistMan/HistMan.h"
00016 #include "MCNtuple/NtpMCRecord.h"
00017 #include "MCNtuple/NtpMCTruth.h"
00018 #include "TruthHelperNtuple/NtpTHRecord.h"
00019 #include "CandNtupleSR/NtpSRRecord.h"
00020 #include "CandNtupleSR/NtpSREvent.h"
00021 #include "CandNtupleSR/NtpSRTrack.h"
00022 #include "CandNtupleSR/NtpSRShower.h"
00023 #include "NueAna/SntpHelpers.h"
00024 
00025 
00026 JOBMODULE(NueReadTJPID, "NueReadTJPID",
00027           "Reads in ana_nue ntuple and pid ntuple");
00028 CVSID("$Id:");
00029 //......................................................................
00030 
00031 NueReadTJPID::NueReadTJPID()
00032 {}
00033 
00034 //......................................................................
00035 
00036 NueReadTJPID::~NueReadTJPID()
00037 {}
00038 
00039 //......................................................................
00040 
00041 JobCResult NueReadTJPID::Ana(const MomNavigator* mom)
00042 {
00043    //get all NueRecords from mom 
00044    //may have more than one per go since mom reads in a snarl's worth of data
00045    //so, this is a little more complicated than just asking for a NueRecord
00046    MSG("NueReadTJPID",Msg::kDebug)<<"***********IN ANA*****************"<<endl;
00047    
00048    TObject *obj=0;
00049    TIter objiter = mom->FragmentIter();
00050    NtpSRRecord *sr=0;
00051    NtpMCRecord *mc=0;
00052    NtpTHRecord *th=0;
00053    vector<NuePID *> vpid;
00054    while((obj=objiter.Next())){
00055       const char *cn=obj->ClassName();
00056 //      MSG("NueReadTJPID",Msg::kDebug)<<"Found a "<<cn<<"."<<endl;
00057       if(strcmp(cn,"NtpSRRecord")==0){
00058          sr = dynamic_cast<NtpSRRecord *>(obj);
00059       }
00060       else if(strcmp(cn,"NtpMCRecord")==0){
00061          mc = dynamic_cast<NtpMCRecord *>(obj);
00062       }
00063       else if(strcmp(cn,"NtpTHRecord")==0){
00064          th = dynamic_cast<NtpTHRecord *>(obj);
00065       }
00066       else if(strcmp(cn,"NuePID")==0){
00067          NuePID *npid  = dynamic_cast<NuePID *>(obj);
00068 //       MSG("NueReadTJPID",Msg::kDebug)<<"Found a NuePID in MOM"
00069 //                                     <<" Snarl "<<npid->GetHeader().GetSnarl()
00070 //                                     <<" Event "<<npid->GetHeader().GetEventNo()<<endl;
00071          vpid.push_back(npid);
00072       }
00073       else{
00074          continue;
00075       }
00076    }
00077    if(!sr){
00078       MSG("NueReadTJPID",Msg::kError)<<"No NtpSRRecord in Mom"<<endl;
00079       return JobCResult::kFailed;
00080    }
00081    
00082    if(!mc){
00083       MSG("NueReadTJPID",Msg::kError)<<"No NtpMCRecord in Mom"<<endl;
00084       return JobCResult::kFailed;
00085    }   
00086    if(!th){
00087       MSG("NueReadTJPID",Msg::kError)<<"No NtpTHRecord in Mom"<<endl;
00088       return JobCResult::kFailed;
00089    }
00090    
00091 //   MSG("NueReadTJPID",Msg::kDebug)<<"SR snarl: "<<sr->GetHeader().GetSnarl()<<endl;
00092 
00093    if(vpid.size()==0){
00094       MSG("NueReadTJPID",Msg::kError)<<"No NuePID Records in Mom"<<endl;
00095       return JobCResult::kFailed;
00096    }
00097       
00098    //use HistMan to plot something for events that pass/fail
00099 
00100    static HistMan *hm = new HistMan("tjpid");
00101    static TH1F *htEall=0;
00102    static TH1F *htEsig=0;
00103    static TH1F *htEnumu=0;
00104    static TH1F *htEnutau=0;
00105    static TH1F *htEbe=0;
00106    static TH1F *htEnc=0;
00107 
00108    static TH1F *hphall=0;
00109    static TH1F *hphsig=0;
00110    static TH1F *hphnumu=0;
00111    static TH1F *hphnutau=0;
00112    static TH1F *hphbe=0;
00113    static TH1F *hphnc=0;
00114 
00115    static TH1F *htrklnall=0;
00116    static TH1F *htrklnsig=0;
00117    static TH1F *htrklnnumu=0;
00118    static TH1F *htrklnnutau=0;
00119    static TH1F *htrklnbe=0;
00120    static TH1F *htrklnnc=0;
00121 
00122    static TH1F *hshwlnall=0;
00123    static TH1F *hshwlnsig=0;
00124    static TH1F *hshwlnnumu=0;
00125    static TH1F *hshwlnnutau=0;
00126    static TH1F *hshwlnbe=0;
00127    static TH1F *hshwlnnc=0;
00128 
00129 
00130 
00131    htEall = hm->Book<TH1F>("htEall","htEall",50,0,50);
00132    htEsig = hm->Book<TH1F>("htEsig","htEsig",50,0,50);
00133    htEnumu = hm->Book<TH1F>("htEnumu","htEnumu",50,0,50);
00134    htEnutau = hm->Book<TH1F>("htEnutau","htEnutau",50,0,50);
00135    htEbe = hm->Book<TH1F>("htEbe","htEbe",50,0,50);
00136    htEnc = hm->Book<TH1F>("htEnc","htEnc",50,0,50);
00137 
00138    hphall = hm->Book<TH1F>("hphall","hphall",50,0,50);
00139    hphsig = hm->Book<TH1F>("hphsig","hphsig",50,0,50);
00140    hphnumu = hm->Book<TH1F>("hphnumu","hphnumu",50,0,50);
00141    hphnutau = hm->Book<TH1F>("hphnutau","hphnutau",50,0,50);
00142    hphbe = hm->Book<TH1F>("hphbe","hphbe",50,0,50);
00143    hphnc = hm->Book<TH1F>("hphnc","hphnc",50,0,50);
00144 
00145    htrklnall = hm->Book<TH1F>("htrklnall","htrklnall",50,0,50);
00146    htrklnsig = hm->Book<TH1F>("htrklnsig","htrklnsig",50,0,50);
00147    htrklnnumu = hm->Book<TH1F>("htrklnnumu","htrklnnumu",50,0,50);
00148    htrklnnutau = hm->Book<TH1F>("htrklnnutau","htrklnnutau",50,0,50);
00149    htrklnbe = hm->Book<TH1F>("htrklnbe","htrklnbe",50,0,50);
00150    htrklnnc = hm->Book<TH1F>("htrklnnc","htrklnnc",50,0,50);
00151 
00152    hshwlnall = hm->Book<TH1F>("hshwlnall","hshwlnall",50,0,50);
00153    hshwlnsig = hm->Book<TH1F>("hshwlnsig","hshwlnsig",50,0,50);
00154    hshwlnnumu = hm->Book<TH1F>("hshwlnnumu","hshwlnnumu",50,0,50);
00155    hshwlnnutau = hm->Book<TH1F>("hshwlnnutau","hshwlnnutau",50,0,50);
00156    hshwlnbe = hm->Book<TH1F>("hshwlnbe","hshwlnbe",50,0,50);
00157    hshwlnnc = hm->Book<TH1F>("hshwlnnc","hshwlnnc",50,0,50);
00158 
00159 
00160    
00161    //so, mom will match up snarls for us,
00162    //but, we have to match up events for ourselves.
00163    for(unsigned int i=0;i<vpid.size();i++){
00164       int evtno = vpid[i]->GetHeader().GetEventNo();
00165       int mcindex=0;
00166       if(evtno<0){
00167          MSG("NueReadTJPID",Msg::kDebug)<<"can not get mctruth for event "<<evtno<<endl;
00168       }
00169       else{
00170          MSG("NueReadTJPID",Msg::kDebug)<<"Trying to get mc index "<<evtno<<endl;
00171          mcindex = SntpHelpers::GetEvent2MCIndex(evtno,th);
00172          MSG("NueReadTJPID",Msg::kDebug)<<"got mc index "<<mcindex<<endl;
00173       }
00174       NtpMCTruth *mcth = SntpHelpers::GetMCTruth(mcindex,mc);
00175       if(mcth==0){
00176          MSG("NueReadTJPID",Msg::kError)<<"can not get mctruth for event "<<evtno<<endl;
00177          continue;
00178       }
00179 
00180       NtpSREvent *event = 0;
00181       if(evtno>=0){
00182          event = SntpHelpers::GetEvent(evtno,sr);
00183       }
00184       //loop over tracks in this event, find longest
00185       int longtrack=0;
00186       if(event!=0){
00187          for(int j=0;j<event->ntrack;j++){
00188             int tindex = SntpHelpers::GetTrackIndex(j,event);
00189             NtpSRTrack *track = SntpHelpers::GetTrack(tindex,sr);
00190             if(longtrack<track->plane.n){
00191                longtrack = track->plane.n;
00192             }
00193          }
00194       }
00195       //loop over showers in this event, find longest
00196       int longshower=0;
00197       if(event!=0){
00198          for(int j=0;j<event->nshower;j++){
00199             int sindex = SntpHelpers::GetShowerIndex(j,event);
00200             NtpSRShower *shower = SntpHelpers::GetShower(sindex,sr);
00201             if(longshower<shower->plane.n){
00202                longshower = shower->plane.n;
00203             }
00204          }
00205       }
00206 
00207       hm->Fill1d("htEall",mcth->p4neu[3]);
00208       if(event!=0){
00209          hm->Fill1d("hphall",event->ph.pe);
00210       }
00211       hm->Fill1d("htrklnall",longtrack);
00212       hm->Fill1d("hshwlnall",longshower);
00213       
00214       int pass = vpid[i]->IsNue;
00215       if(pass==1){
00216          int cls=0;
00217          if(mcth->iaction==0){
00218             cls=5;
00219          }
00220          else if(abs(mcth->inu)==12){
00221             if(abs(mcth->inunoosc)==12){
00222                cls=4;
00223             }
00224             else if(abs(mcth->inu)==14){
00225                cls=1;
00226             }
00227          }
00228          else if(abs(mcth->inu)==14){
00229             cls=2;
00230          }
00231          else if(abs(mcth->inu)==16){
00232             cls=3;
00233          }
00234          if(cls==1){
00235             hm->Fill1d("htEsig",mcth->p4neu[3]);
00236             if(event!=0) hm->Fill1d("hphsig",event->ph.pe);
00237             hm->Fill1d("htrklnsig",longtrack);
00238             hm->Fill1d("hshwlnsig",longshower);
00239          }
00240          else if(cls==2){
00241             hm->Fill1d("htEnumu",mcth->p4neu[3]);
00242             if(event!=0) hm->Fill1d("hphnumu",event->ph.pe);
00243             hm->Fill1d("htrklnnumu",longtrack);
00244             hm->Fill1d("hshwlnnumu",longshower);
00245          }
00246          else if(cls==3){
00247             hm->Fill1d("htEnutau",mcth->p4neu[3]);
00248             if(event!=0) hm->Fill1d("hphnutau",event->ph.pe);
00249             hm->Fill1d("htrklnnutau",longtrack);
00250             hm->Fill1d("hshwlnnutau",longshower);
00251          }
00252          else if(cls==4){
00253             hm->Fill1d("htEbe",mcth->p4neu[3]);
00254             if(event!=0) hm->Fill1d("hphbe",event->ph.pe);
00255             hm->Fill1d("htrklnbe",longtrack);
00256             hm->Fill1d("hshwlnbe",longshower);
00257          }
00258          else if(cls==5){
00259             hm->Fill1d("htEnc",mcth->p4neu[3]);
00260             if(event!=0) hm->Fill1d("hphnc",event->ph.pe);
00261             hm->Fill1d("htrklnnc",longtrack);
00262             hm->Fill1d("hshwlnnc",longshower);
00263          }
00264       }
00265    }
00266 
00267    return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00268 }
00269