loon: CompareMST.cxx Source File

00001 
00002 // $Id: CompareMST.cxx,v 1.2 2008/11/19 18:22:51 rhatcher 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 "Conventions/Detector.h"
00011 #include "NueAna/CompareMST.h"
00012 #include "NueAna/NueRecord.h"
00013 #include "MessageService/MsgService.h"
00014 #include "MinosObjectMap/MomNavigator.h"
00015 #include "MessageService/MsgService.h"
00016 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00017 #include "HistMan/HistMan.h"
00018 
00019 const int CALend=110;
00020 const int SM1end=230;
00021 const int SM2end=484;
00022 
00023 JOBMODULE(CompareMST, "CompareMST",
00024           "Makes hists to compare MST variables between detectors and truth classes");
00025 CVSID("$Id:");
00026 //......................................................................
00027 
00028 CompareMST::CompareMST():
00029   counter(0)
00030 {}
00031 
00032 //......................................................................
00033 
00034 CompareMST::~CompareMST()
00035 {}
00036 
00037 //......................................................................
00038 void CompareMST::BeginJob()
00039 {
00040 
00041   if(counter%1000==0){
00042     cout<<"On entry "<<counter<<endl;
00043   }
00044   counter++;
00045  
00046  static HistMan *hm = new HistMan("mstcomp");
00047   for(int d=0;d<2;d++){
00048     char det[2];
00049     if(d==0){
00050           sprintf(det,"f");
00051     }
00052     else if(d==1){
00053       sprintf(det,"n");
00054     }
00055     else{
00056       sprintf(det,"u");
00057     }
00058     for(int b=0;b<4;b++){
00059       char bgc[6];
00060       if(b==0){
00061         sprintf(bgc,"nue");
00062       }
00063       else if(b==1){
00064         sprintf(bgc,"numu");
00065       }
00066       else if(b==2){
00067         sprintf(bgc,"nutau");
00068       }
00069       else if(b==3){
00070         sprintf(bgc,"nc");
00071       }
00072       else{
00073         sprintf(bgc,"u");
00074       }
00075       for(int r=0;r<4;r++){
00076         char id[20];
00077         sprintf(id,"%s_%s_%d",det,bgc,r+1001);
00078         char nemsg[50];
00079         char newtot[50];
00080         char nentot[50];
00081         char new1[50];
00082         char nen1[50];
00083         char new1nn1[50];
00084         char nesmtot[50];
00085         char nesm1[50];
00086         char neal[50];
00087         char nep[50];
00088         char new4[50];
00089         char nen4[50];
00090         char new4n4[50];
00091         char nesm4w4[50];
00092         char newall[50];
00093         char nemall[50];
00094         char neb1[50];
00095         char neb25[50];
00096         char nebranch[50];
00097 
00098         sprintf(nemsg,"enmsg_%s",id);
00099         sprintf(newtot,"ewtot_%s",id);
00100         sprintf(nentot,"entot_%s",id);
00101         sprintf(new1,"ew1_%s",id);
00102         sprintf(nen1,"en1_%s",id);
00103         sprintf(new1nn1,"ew1nn1_%s",id);
00104         sprintf(nesmtot,"esmtot_%s",id);
00105         sprintf(nesm1,"esm1_%s",id);
00106         sprintf(neal,"eal_%s",id);
00107         sprintf(nep,"ep_%s",id);
00108         sprintf(new4,"ew4_%s",id);
00109         sprintf(nen4,"en4_%s",id);
00110         sprintf(new4n4,"ew4n4_%s",id);
00111         sprintf(nesm4w4,"esm4w4_%s",id);
00112         sprintf(newall,"ewall_%s",id);
00113         sprintf(nemall,"emall_%s",id);
00114         sprintf(neb1,"eb1_%s",id);
00115         sprintf(neb25,"eb25_%s",id);
00116         sprintf(nebranch,"enbranch_%s",id);
00117 
00118         /*
00119         TH1F *hemsg=hm->Book<TH1F>(nemsg,nemsg,11,-0.5,10.5);
00120         TH1F *hewtot=hm->Book<TH1F>(newtot,newtot,100,0,2000);
00121         TH1F *hentot=hm->Book<TH1F>(nentot,nentot,100,0,300);
00122         TH1F *hew1=hm->Book<TH1F>(new1,new1,100,0,2000);
00123         TH1F *hen1=hm->Book<TH1F>(nen1,nen1,100,0,300);
00124         TH1F *hew1nn1=hm->Book<TH1F>(new1nn1,new1nn1,100,0,50);
00125         TH1F *hesmtot=hm->Book<TH1F>(nesmtot,nesmtot,100,0,5000);
00126         TH1F *hesm1=hm->Book<TH1F>(nesm1,nesm1,100,0,5000);
00127         TH1F *heal=hm->Book<TH1F>(neal,neal,51,-0.01,1.01);
00128         TH1F *hep=hm->Book<TH1F>(nep,nep,50,0,100);
00129         TH1F *hew4=hm->Book<TH1F>(new4,new4,50,0,1000);
00130         TH1F *hen4=hm->Book<TH1F>(nen4,nen4,50,0,100);
00131         TH1F *hew4n4=hm->Book<TH1F>(new4n4,new4n4,100,0,50);
00132         TH1F *hesm4w4=hm->Book<TH1F>(nesm4w4,nesm4w4,100,0,5);
00133         TH1F *hewall=hm->Book<TH1F>(newall,newall,100,0,100);
00134         TH1F *hemall=hm->Book<TH1F>(nemall,nemall,100,0,100);
00135         TH1F *heb1=hm->Book<TH1F>(neb1,neb1,10,0,1);
00136         TH1F *heb25=hm->Book<TH1F>(neb25,neb25,10,0,1);
00137         TH1F *hebranch=hm->Book<TH1F>(nebranch,nebranch,21,-0.5,20.5);
00138         */
00139         hm->Book<TH1F>(nemsg,nemsg,11,-0.5,10.5);
00140         hm->Book<TH1F>(newtot,newtot,100,0,2000);
00141         hm->Book<TH1F>(nentot,nentot,100,0,300);
00142         hm->Book<TH1F>(new1,new1,100,0,2000);
00143         hm->Book<TH1F>(nen1,nen1,100,0,300);
00144         hm->Book<TH1F>(new1nn1,new1nn1,100,0,50);
00145         hm->Book<TH1F>(nesmtot,nesmtot,100,0,5000);
00146         hm->Book<TH1F>(nesm1,nesm1,100,0,5000);
00147         hm->Book<TH1F>(neal,neal,51,-0.01,1.01);
00148         hm->Book<TH1F>(nep,nep,50,0,100);
00149         hm->Book<TH1F>(new4,new4,50,0,1000);
00150         hm->Book<TH1F>(nen4,nen4,50,0,100);
00151         hm->Book<TH1F>(new4n4,new4n4,100,0,50);
00152         hm->Book<TH1F>(nesm4w4,nesm4w4,100,0,5);
00153         hm->Book<TH1F>(newall,newall,100,0,100);
00154         hm->Book<TH1F>(nemall,nemall,100,0,100);
00155         hm->Book<TH1F>(neb1,neb1,10,0,1);
00156         hm->Book<TH1F>(neb25,neb25,10,0,1);
00157         hm->Book<TH1F>(nebranch,nebranch,21,-0.5,20.5);
00158 
00159 
00160         char nomsg[50];
00161         char nowtot[50];
00162         char nontot[50];
00163         char now1[50];
00164         char non1[50];
00165         char now1nn1[50];
00166         char nosmtot[50];
00167         char nosm1[50];
00168         char noal[50];
00169         char nop[50];
00170         char now4[50];
00171         char non4[50];
00172         char now4n4[50];
00173         char nosm4w4[50];
00174         char nowall[50];
00175         char nomall[50];
00176         char nob1[50];
00177         char nob25[50];
00178         char nobranch[50];
00179 
00180         sprintf(nomsg,"onmsg_%s",id);
00181         sprintf(nowtot,"owtot_%s",id);
00182         sprintf(nontot,"ontot_%s",id);
00183         sprintf(now1,"ow1_%s",id);
00184         sprintf(non1,"on1_%s",id);
00185         sprintf(now1nn1,"ow1nn1_%s",id);
00186         sprintf(nosmtot,"osmtot_%s",id);
00187         sprintf(nosm1,"osm1_%s",id);
00188         sprintf(noal,"oal_%s",id);
00189         sprintf(nop,"op_%s",id);
00190         sprintf(now4,"ow4_%s",id);
00191         sprintf(non4,"on4_%s",id);
00192         sprintf(now4n4,"ow4n4_%s",id);
00193         sprintf(nosm4w4,"osm4w4_%s",id);
00194         sprintf(nowall,"owall_%s",id);
00195         sprintf(nomall,"omall_%s",id);
00196         sprintf(nob1,"ob1_%s",id);
00197         sprintf(nob25,"ob25_%s",id);
00198         sprintf(nobranch,"onbranch_%s",id);
00199         /*
00200         TH1F *homsg=hm->Book<TH1F>(nomsg,nomsg,11,-0.5,10.5);
00201         TH1F *howtot=hm->Book<TH1F>(nowtot,nowtot,100,0,2000);
00202         TH1F *hontot=hm->Book<TH1F>(nontot,nontot,100,0,300);
00203         TH1F *how1=hm->Book<TH1F>(now1,now1,100,0,2000);
00204         TH1F *hon1=hm->Book<TH1F>(non1,non1,100,0,300);
00205         TH1F *how1nn1=hm->Book<TH1F>(now1nn1,now1nn1,100,0,50);
00206         TH1F *hosmtot=hm->Book<TH1F>(nosmtot,nosmtot,100,0,5000);
00207         TH1F *hosm1=hm->Book<TH1F>(nosm1,nosm1,100,0,5000);
00208         TH1F *hoal=hm->Book<TH1F>(noal,noal,51,-0.01,1.01);
00209         TH1F *hop=hm->Book<TH1F>(nop,nop,50,0,100);
00210         TH1F *how4=hm->Book<TH1F>(now4,now4,50,0,1000);
00211         TH1F *hon4=hm->Book<TH1F>(non4,non4,50,0,100);
00212         TH1F *how4n4=hm->Book<TH1F>(now4n4,now4n4,100,0,50);
00213         TH1F *hosm4w4=hm->Book<TH1F>(nosm4w4,nosm4w4,100,0,5);
00214         TH1F *howall=hm->Book<TH1F>(nowall,nowall,100,0,100);
00215         TH1F *homall=hm->Book<TH1F>(nomall,nomall,100,0,100);
00216         TH1F *hob1=hm->Book<TH1F>(nob1,nob1,10,0,1);
00217         TH1F *hob25=hm->Book<TH1F>(nob25,nob25,10,0,1);
00218         TH1F *hobranch=hm->Book<TH1F>(nobranch,nobranch,21,-0.5,20.5);
00219         */      
00220         hm->Book<TH1F>(nomsg,nomsg,11,-0.5,10.5);
00221         hm->Book<TH1F>(nowtot,nowtot,100,0,2000);
00222         hm->Book<TH1F>(nontot,nontot,100,0,300);
00223         hm->Book<TH1F>(now1,now1,100,0,2000);
00224         hm->Book<TH1F>(non1,non1,100,0,300);
00225         hm->Book<TH1F>(now1nn1,now1nn1,100,0,50);
00226         hm->Book<TH1F>(nosmtot,nosmtot,100,0,5000);
00227         hm->Book<TH1F>(nosm1,nosm1,100,0,5000);
00228         hm->Book<TH1F>(noal,noal,51,-0.01,1.01);
00229         hm->Book<TH1F>(nop,nop,50,0,100);
00230         hm->Book<TH1F>(now4,now4,50,0,1000);
00231         hm->Book<TH1F>(non4,non4,50,0,100);
00232         hm->Book<TH1F>(now4n4,now4n4,100,0,50);
00233         hm->Book<TH1F>(nosm4w4,nosm4w4,100,0,5);
00234         hm->Book<TH1F>(nowall,nowall,100,0,100);
00235         hm->Book<TH1F>(nomall,nomall,100,0,100);
00236         hm->Book<TH1F>(nob1,nob1,10,0,1);
00237         hm->Book<TH1F>(nob25,nob25,10,0,1);
00238         hm->Book<TH1F>(nobranch,nobranch,21,-0.5,20.5);
00239 
00240       }
00241     }
00242   }
00243 }
00244 
00245 
00246 JobCResult CompareMST::Ana(const MomNavigator* mom)
00247 {
00248    //get all NueRecords from mom 
00249    //may have more than one per go since mom reads in a snarl's worth of data
00250    //so, this is a little more complicated than just asking for a NueRecord
00251    TObject *obj=0;
00252    TIter objiter = mom->FragmentIter();
00253    while((obj=objiter.Next())){
00254       NueRecord *nr = static_cast<NueRecord *>(obj);
00255       if(nr){
00256          MSG("CompareMST",Msg::kDebug)<<"Found a NueRecord in MOM"<<endl;
00257       }
00258       else{
00259          MSG("CompareMST",Msg::kError)<<"Didn't find a NueRecord in MOM"<<endl;
00260          continue;
00261       }
00262       MSG("CompareMST",Msg::kDebug)<<"Found a NueRecord "<<nr<<endl;
00263       
00264       //need at least 1 sucessfully reconstructed event
00265       if(nr->GetHeader().GetEventNo()<0){
00266         continue;
00267       }
00268       //only look at events for which mst actually gets calculated
00269       if(nr->GetHeader().GetTrackLength()>25){
00270         continue;
00271       }
00272       //only look at contained events (vertex contained, end doesn't go past 
00273       //calorimeter or sm gap
00274       if(nr->GetHeader().GetVldContext().GetDetector()==Detector::kFar){
00275         if(sqrt(pow(nr->srevent.vtxX,2)+pow(nr->srevent.vtxY,2))>3.87){
00276           continue;
00277         }
00278         if(nr->srevent.vtxZ<1||(nr->srevent.vtxZ>14&&nr->srevent.vtxZ<17)||
00279            nr->srevent.vtxZ>20){
00280           continue;
00281         }
00282         if((nr->srevent.begPlane<SM1end&&nr->srevent.endPlane>SM1end-2)||
00283            (nr->srevent.begPlane>SM1end&&nr->srevent.endPlane>SM2end-2)){
00284           continue;
00285         }         
00286       }
00287       else if(nr->GetHeader().GetVldContext().GetDetector()==Detector::kNear){
00288         if(sqrt(pow(nr->srevent.vtxX-1.5,2)+pow(nr->srevent.vtxY,2))>1.0){
00289           continue;
00290         }
00291         if(nr->srevent.vtxZ<4||nr->srevent.vtxZ>6.5){
00292           continue;
00293         }
00294         if(nr->srevent.endPlane>CALend-2){
00295           continue;
00296         }         
00297       }
00298       //only look at events where mst output makes sense, nn1>3
00299       if(nr->mstvars.enn1<=3||nr->mstvars.onn1<=3){
00300         continue;
00301       }
00302       
00303       //make id string
00304       Detector::Detector_t d = nr->GetHeader().GetVldContext().GetDetector();
00305       char evdet[2];
00306       if(d==Detector::kFar){
00307         sprintf(evdet,"f");
00308       }
00309       else if(d==Detector::kNear){
00310         sprintf(evdet,"n");
00311       }
00312       else{
00313         sprintf(evdet,"u");
00314       }
00315       
00316       char evbgc[6];
00317       if(nr->mctrue.interactionType==1){
00318         if(abs(nr->mctrue.nuFlavor)==12){
00319           sprintf(evbgc,"nue");
00320         }
00321         else if(abs(nr->mctrue.nuFlavor)==14){
00322           sprintf(evbgc,"numu");
00323         }
00324         else if(abs(nr->mctrue.nuFlavor)==16){
00325           sprintf(evbgc,"nutau");
00326         }
00327         else{
00328           sprintf(evbgc,"u");
00329         }
00330       }
00331       else{
00332         sprintf(evbgc,"nc");
00333       }
00334 
00335       char evid[20];
00336       sprintf(evid,"%s_%s_%d",evdet,evbgc,nr->mctrue.resonanceCode);
00337 
00338       //make histname strings
00339       char nemsg[50];
00340       char newtot[50];
00341       char nentot[50];
00342       char new1[50];
00343       char nen1[50];
00344       char new1nn1[50];
00345       char nesmtot[50];
00346       char nesm1[50];
00347       char neal[50];
00348       char nep[50];
00349       char new4[50];
00350       char nen4[50];
00351       char new4n4[50];
00352       char nesm4w4[50];
00353       char newall[50];
00354       char nemall[50];
00355       char neb1[50];
00356       char neb25[50];
00357       char nebranch[50];
00358       
00359       sprintf(nemsg,"enmsg_%s",evid);
00360       sprintf(newtot,"ewtot_%s",evid);
00361       sprintf(nentot,"entot_%s",evid);
00362       sprintf(new1,"ew1_%s",evid);
00363       sprintf(nen1,"en1_%s",evid);
00364       sprintf(new1nn1,"ew1nn1_%s",evid);
00365       sprintf(nesmtot,"esmtot_%s",evid);
00366       sprintf(nesm1,"esm1_%s",evid);
00367       sprintf(neal,"eal_%s",evid);
00368       sprintf(nep,"ep_%s",evid);
00369       sprintf(new4,"ew4_%s",evid);
00370       sprintf(nen4,"en4_%s",evid);
00371       sprintf(new4n4,"ew4n4_%s",evid);
00372       sprintf(nesm4w4,"esm4w4_%s",evid);
00373       sprintf(newall,"ewall_%s",evid);
00374       sprintf(nemall,"emall_%s",evid);
00375       sprintf(neb1,"eb1_%s",evid);
00376       sprintf(neb25,"eb25_%s",evid);
00377       sprintf(nebranch,"enbranch_%s",evid);
00378 
00379       char nomsg[50];
00380       char nowtot[50];
00381       char nontot[50];
00382       char now1[50];
00383       char non1[50];
00384       char now1nn1[50];
00385       char nosmtot[50];
00386       char nosm1[50];
00387       char noal[50];
00388       char nop[50];
00389       char now4[50];
00390       char non4[50];
00391       char now4n4[50];
00392       char nosm4w4[50];
00393       char nowall[50];
00394       char nomall[50];
00395       char nob1[50];
00396       char nob25[50];
00397       char nobranch[50];
00398 
00399       sprintf(nomsg,"onmsg_%s",evid);
00400       sprintf(nowtot,"owtot_%s",evid);
00401       sprintf(nontot,"ontot_%s",evid);
00402       sprintf(now1,"ow1_%s",evid);
00403       sprintf(non1,"on1_%s",evid);
00404       sprintf(now1nn1,"ow1nn1_%s",evid);
00405       sprintf(nosmtot,"osmtot_%s",evid);
00406       sprintf(nosm1,"osm1_%s",evid);
00407       sprintf(noal,"oal_%s",evid);
00408       sprintf(nop,"op_%s",evid);
00409       sprintf(now4,"ow4_%s",evid);
00410       sprintf(non4,"on4_%s",evid);
00411       sprintf(now4n4,"ow4n4_%s",evid);
00412       sprintf(nosm4w4,"osm4w4_%s",evid);
00413       sprintf(nowall,"owall_%s",evid);
00414       sprintf(nomall,"omall_%s",evid);
00415       sprintf(nob1,"ob1_%s",evid);
00416       sprintf(nob25,"ob25_%s",evid);
00417       sprintf(nobranch,"onbranch_%s",evid);
00418 
00419       //fill histograms      
00420       static HistMan *hm = new HistMan("mstcomp");
00421       //      cout<<"Printing hm "<<endl;
00422       //      hm->BaseFolder().ls();
00423       //      cout<<"******************************"<<endl;
00424 
00425       hm->Fill1d(nemsg,nr->mstvars.enmsg);
00426       hm->Fill1d(newtot,nr->mstvars.ewtot);
00427       hm->Fill1d(nentot,nr->mstvars.enntot);
00428       hm->Fill1d(new1,nr->mstvars.ew1);
00429       hm->Fill1d(nen1,nr->mstvars.enn1);
00430       if(nr->mstvars.enn1!=0){
00431         hm->Fill1d(new1nn1,nr->mstvars.ew1/nr->mstvars.enn1);
00432       }
00433       hm->Fill1d(nesmtot,nr->mstvars.esmtot);
00434       hm->Fill1d(nesm1,nr->mstvars.esm1);
00435       hm->Fill1d(neal,nr->mstvars.ealpha);
00436       hm->Fill1d(nep,nr->mstvars.eeprob);
00437       hm->Fill1d(new4,nr->mstvars.e4w);
00438       hm->Fill1d(nen4,nr->mstvars.e4nn);
00439       if(nr->mstvars.e4nn!=0){
00440         hm->Fill1d(new4n4,nr->mstvars.e4w/nr->mstvars.e4nn);
00441       }
00442       if(nr->mstvars.e4w!=0){
00443         hm->Fill1d(nesm4w4,nr->mstvars.e4sm/nr->mstvars.e4w);
00444       }
00445       for(int i=0;i<nr->mstvars.enn1;i++){
00446         hm->Fill1d(newall,nr->mstvars.eallw1[i]);
00447         hm->Fill1d(nemall,nr->mstvars.eallm1[i]);
00448       }
00449       hm->Fill1d(neb1,nr->mstvars.eb1);
00450       hm->Fill1d(neb25,nr->mstvars.eb25);
00451       hm->Fill1d(nebranch,nr->mstvars.enbranch);
00452 
00453       hm->Fill1d(nomsg,nr->mstvars.onmsg);
00454       hm->Fill1d(nowtot,nr->mstvars.owtot);
00455       hm->Fill1d(nontot,nr->mstvars.onntot);
00456       hm->Fill1d(now1,nr->mstvars.ow1);
00457       hm->Fill1d(non1,nr->mstvars.onn1);
00458       if(nr->mstvars.onn1!=0){
00459         hm->Fill1d(now1nn1,nr->mstvars.ow1/nr->mstvars.onn1);
00460       }
00461       hm->Fill1d(nosmtot,nr->mstvars.osmtot);
00462       hm->Fill1d(nosm1,nr->mstvars.osm1);
00463       hm->Fill1d(noal,nr->mstvars.oalpha);
00464       hm->Fill1d(nop,nr->mstvars.oeprob);
00465       hm->Fill1d(now4,nr->mstvars.o4w);
00466       hm->Fill1d(non4,nr->mstvars.o4nn);
00467       if(nr->mstvars.o4nn!=0){
00468         hm->Fill1d(now4n4,nr->mstvars.o4w/nr->mstvars.o4nn);
00469       }
00470       if(nr->mstvars.o4w!=0){
00471         hm->Fill1d(nosm4w4,nr->mstvars.o4sm/nr->mstvars.o4w);
00472       }
00473       for(int i=0;i<nr->mstvars.onn1;i++){
00474         hm->Fill1d(nowall,nr->mstvars.oallw1[i]);
00475         hm->Fill1d(nomall,nr->mstvars.oallm1[i]);
00476       }
00477       hm->Fill1d(nob1,nr->mstvars.ob1);
00478       hm->Fill1d(nob25,nr->mstvars.ob25);
00479       hm->Fill1d(nobranch,nr->mstvars.onbranch);
00480 
00481 
00482 
00483    }
00484    return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00485 }
00486 
00488 void CompareMST::EndJob()
00489 {}