loon: AnalysisInfoAna.cxx Source File

00001 
00016 #include <cassert>
00017 
00018 #include "TSystem.h"
00019 #include "CandNtupleSR/NtpSRRecord.h"
00020 #include "CandNtupleSR/NtpSREvent.h"
00021 #include "MessageService/MsgService.h"
00022 #include "StandardNtuple/NtpStRecord.h"
00023 #include "NueAna/NueAnaTools/SntpHelpers.h"
00024 #include "NueAna/AnalysisInfoAna.h"
00025 #include "NueAna/NueAnalysisCuts.h"
00026 #include "AnalysisNtuples/ANtpDefaultValue.h"
00027 #include "Mad/MadNsID.h"
00028 #include "Mad/MadDpID.h"
00029 #include "Mad/MadAbID.h"
00030 #include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
00031 #include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
00032 #include "DataUtil/infid.h"
00033 
00034 #include "NueAna/Module/SetKNNModule.h"  
00035 #include "PhysicsNtuple/Handle.h"
00036 #include "PhysicsNtuple/Factory.h"
00037   
00038 
00039                                                                                 
00040 MadNsID AnalysisInfoAna::nsid;
00041 MadDpID AnalysisInfoAna::dpid;
00042 MadAbID AnalysisInfoAna::abid;
00043 
00044 bool AnalysisInfoAna::readabidfile=false;
00045                                                                                 
00046 CVSID("$Id: AnalysisInfoAna.cxx,v 1.24 2009/10/13 16:06:03 scavan Exp $");
00047 
00048 AnalysisInfoAna::AnalysisInfoAna(AnalysisInfoNue &anai):
00049    fAnalysisInfo(anai)
00050 {
00051     fDetectorType =  Detector::kUnknown;
00052 }
00053                                                                                 
00054 AnalysisInfoAna::~AnalysisInfoAna()
00055 {}
00056                                                                                 
00057 void AnalysisInfoAna::Analyze(int evtn, RecRecordImp<RecCandHeader> *srobj)
00058 {
00059   MSG("AnalysisInfoAna",Msg::kDebug)
00060       <<"Entering AnalysisInfoAna::Analyze"<<endl;
00061                                                                                 
00062   if(srobj==0){
00063       return;
00064   }  if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00065         ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00066      return;
00067   }
00068                                                                                 
00069   const RecCandHeader *ntpHeader = &(srobj->GetHeader());
00070   fDetectorType = ntpHeader->GetVldContext().GetDetector();
00071   
00072   //Not applicable to CalDet:
00073   if(fDetectorType==Detector::kCalDet) return;
00074   
00075   bool foundst=false;
00076 
00077   ANtpRecoNtpManipulator *ntpManipulator = 0;
00078                                                                                 
00079   NtpSRRecord *sr = 0;
00080   NtpStRecord *st=dynamic_cast<NtpStRecord *>(srobj);
00081 
00082   NtpSREvent *event = 0;
00083   NtpSRTrack *track = 0;
00084   NtpSRShower *shower = 0;
00085   NtpSREventSummary *eventSummary;
00086   SimFlag::SimFlag_t sim;
00087 
00088   if(st != 0)
00089   {
00090       foundst = true;
00091       ntpManipulator =  new ANtpRecoNtpManipulator(st);
00092       eventSummary = &(st->evthdr);
00093       sim = st->GetHeader().GetVldContext().GetSimFlag();
00094   }
00095   else{
00096       sr=dynamic_cast<NtpSRRecord *>(srobj);
00097       ntpManipulator =  new ANtpRecoNtpManipulator(sr);
00098       eventSummary = &(sr->evthdr);
00099       sim = sr->GetHeader().GetVldContext().GetSimFlag();
00100   }
00101 
00102   ntpManipulator->SetPrimaryShowerCriteria(0,1); //nplanes, total ph
00103   ntpManipulator->SetPrimaryTrackCriteria(0,1,0); //nplanes, length, total ph
00104   ANtpEventManipulator * ntpEventManipulator =
00105                             ntpManipulator->GetEventManipulator();
00106                                                                                 
00107   ntpEventManipulator->SetEventInSnarl(evtn);
00108   event  = ntpEventManipulator->GetEvent();
00109   shower = SntpHelpers::GetPrimaryShower(evtn, st);
00110   track  = SntpHelpers::GetPrimaryTrack(evtn,st);
00111 
00112   Float_t vertexX = event->vtx.x; //Munits::meters
00113   Float_t vertexY = event->vtx.y; //Munits::meters
00114   Float_t vertexZ = event->vtx.z; //Munits::meters
00115                                                                                
00116   if(ReleaseType::IsCedar(release)){
00117     NtpVtxFinder vtxf;
00118     vtxf.SetTargetEvent(st, event, track);
00119     if(vtxf.FindVertex() > 0){
00120        vertexX = vtxf.VtxX();
00121        vertexY = vtxf.VtxY();
00122        vertexZ = vtxf.VtxZ();
00123     }
00124   }
00125 
00126                     
00127   fAnalysisInfo.isFullyContained =  
00128      IsFidAll(vertexX, vertexY, vertexZ, event);
00129 
00130   Int_t test = 0;
00131   bool isMC = (sim == SimFlag::kMC);
00132   if(fDetectorType == Detector::kNear)
00133       test = NueConvention::IsInsideNearFiducial_Nue_Standard(vertexX, vertexY, vertexZ, isMC);
00134   if(fDetectorType == Detector::kFar)
00135       test = NueConvention::IsInsideFarFiducial_Nue_Standard(vertexX, vertexY, vertexZ, isMC);
00136   
00137   fAnalysisInfo.inFiducialVolume = test;
00138 
00139 
00140   if(fDetectorType==Detector::kFar ){
00141        dpid.SetPHCorrection(1.018);
00142   }
00143 
00144   BeamType::BeamType_t current_beam = beam;
00145   //  BeamType::BeamType_t current_beam = BeamType::kL250z200i;
00146 
00147   string reco_version = ReleaseType::AsString(release); 
00148   string mc_version = ""; 
00149   if(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC) 
00150     mc_version = ReleaseType::AsString(release); 
00151   if(ReleaseType::IsCarrot(release)) reco_version = "birch"; 
00152   else if(ReleaseType::IsCedar(release)) reco_version = "cedar"; 
00153   else if(ReleaseType::IsDogwood(release)) reco_version = "cedar";  //<-- ! no dogwood file! 
00154   else reco_version = "birch"; 
00155 
00156 
00157   if(dpid.ChoosePDFs(fDetectorType,current_beam,
00158                      reco_version,mc_version))
00159     fAnalysisInfo.dpCCPID = 
00160           dpid.CalcPID(track, event, eventSummary, fDetectorType, 0);
00161 
00162    
00163   if(nsid.ChooseWeightFile(fDetectorType,current_beam)){
00164     if(!nsid.GetPID(event,track,shower,st,fDetectorType,fAnalysisInfo.nsCCPID))
00165        fAnalysisInfo.nsCCPID=ANtpDefVal::kFloat;
00166   }
00167 
00168   if(!readabidfile){
00169     string fname = AnalysisInfoAna::BuildABPIDFile();
00170     abid.ReadPDFs(fname.c_str());    
00171     MSG("AnalysisInfoAna",Msg::kInfo)<<"Reading AB pdfs from "<<fname<<endl;
00172     readabidfile=true;
00173   }
00174 
00175   //now that everything's all set up, compute andy's pid
00176   NtpStRecord *strecord = dynamic_cast<NtpStRecord *>(srobj);
00177   if(strecord==0){}
00178   else{    
00179     fAnalysisInfo.abCCPID=abid.CalcPID(event,strecord);
00180     fAnalysisInfo.abcostheta=abid.GetRecoCosTheta(event,strecord);
00181     fAnalysisInfo.abeventenergy=abid.GetRecoE(event,strecord);
00182     fAnalysisInfo.abtrackcharge=abid.GetTrackEMCharge(event,strecord);
00183     fAnalysisInfo.abtrackenergy=abid.GetTrackPlanes(event,strecord);
00184     fAnalysisInfo.abtracklikeplanes=abid.GetTrackLikePlanes(event,strecord);
00185     fAnalysisInfo.abtrackphfrac=abid.GetTrackPHfrac(event,strecord);
00186     fAnalysisInfo.abtrackphmean=abid.GetTrackPHmean(event,strecord);
00187     fAnalysisInfo.abtrackqpsigmaqp=abid.GetTrackQPsigmaQP(event,strecord);
00188     fAnalysisInfo.aby=abid.GetRecoY(event,strecord);
00189   }
00190 
00191   //rustem's pid
00192    //Loading up Rustems variables
00193    Anp::Handle<StorekNNData> data = Anp::Factory<StorekNNData>::Instance().Get("kNNData");
00194    if(!data.valid())
00195    {
00196       MAXMSG("NueModule", Msg::kError, 1)
00197           << "NueModule::Reco - Handle<StorekNNData> is invalid, assuming no Rustem variable to run"
00198           << endl;
00199    }else{
00200 
00201      float numubar, rel_ang, knn_pid, knn_01, knn_10, knn_20, knn_40; 
00202 
00203      data -> SetPrefix("SNTP");
00204      data -> Get(evtn, "numubar", numubar);
00205      data -> Get(evtn, "rel_ang", rel_ang);
00206      data -> Get(evtn, "knn_pid", knn_pid);
00207      data -> Get(evtn, "knn_01", knn_01);
00208      data -> Get(evtn, "knn_10", knn_10);
00209      data -> Get(evtn, "knn_20", knn_20);
00210      data -> Get(evtn, "knn_40", knn_40);
00211 
00212      fAnalysisInfo.roNScientPlanes = knn_01;
00213      fAnalysisInfo.roMeanTrkSig = knn_10;
00214      fAnalysisInfo.roMeanRatio = knn_20;
00215      fAnalysisInfo.roTrkMeanVsWindow = knn_40;
00216      fAnalysisInfo.roNuMuBar = numubar;
00217      fAnalysisInfo.roRelAng = rel_ang;
00218      fAnalysisInfo.roCCPID  = knn_pid;
00219   }
00220  
00221 
00222   if(ntpManipulator){
00223     delete ntpManipulator;
00224     ntpManipulator=0;
00225   }
00226     
00227 }
00228 
00229 void AnalysisInfoAna::Set3DHit(DeqFloat_t &x
00230                              , DeqFloat_t &y
00231                              , DeqFloat_t &z
00232                              , DeqFloat_t &e){
00233                                                                                 
00234      if(x.size()!=0 && y.size()!=0 && z.size()!=0 && e.size()!=0){
00235                                                                                 
00236            fX=x;
00237            fY=y;
00238            fZ=z;
00239            fE=e;
00240      }
00241 }
00242 
00243 Int_t AnalysisInfoAna::IsFidAll(Float_t vtxX, Float_t vtxY, Float_t vtxZ, NtpSREvent *event)
00244 {
00245   // with 3d hits shut off this function won't do anythnig interesting
00246                                                                                 
00247   Bool_t contained = true;
00248   Int_t test = 0;
00249   if(fX.size()==0 || fY.size()==0 || fZ.size()==0 || fE.size()==0){
00250         return ANtpDefVal::kInt;
00251   }
00252                                                                                 
00253   for(UInt_t i = 0; i < fX.size() && contained; i++){
00254      test = 0;
00255      Float_t x = fX[i] + vtxX;
00256      Float_t y = fY[i] + vtxY;
00257      Float_t z = fZ[i] + vtxZ;
00258      if(fDetectorType == Detector::kNear)
00259       test = NueConvention::IsInsideNearFiducial_Nue_Extended(x,y,z);
00260      if(fDetectorType == Detector::kFar)
00261       test = NueConvention::IsInsideFarFiducial_Nue_Extended(x,y,z);      
00262      if(test <= 0) contained = false;
00263   }
00264                                                                                 
00265   if(event)
00266   if(contained && event->plane.n > 66) //this is approximately 4 meteres which
00267 // is the 3D hit cutoff
00268    {
00269       Float_t x = event->end.x;
00270       Float_t y = event->end.y;
00271       Float_t z = event->end.z;
00272 
00273      if(fDetectorType == Detector::kNear)
00274       test = NueConvention::IsInsideNearFiducial_Nue_Extended(x,y,z);
00275      if(fDetectorType == Detector::kFar)
00276       test = NueConvention::IsInsideFarFiducial_Nue_Extended(x,y,z);
00277      if(test <= 0) contained = false;
00278    }
00279                                                                                 
00280   if(contained) test = 1;
00281   return test;
00282 }
00283 
00284 string AnalysisInfoAna::BuildABPIDFile()
00285 {
00286     string base="";
00287     string pub = getenv("SRT_PUBLIC_CONTEXT");
00288     string priv = getenv("SRT_PRIVATE_CONTEXT");
00289 
00290     if(priv!=""&&priv!="."){
00291       // check if directory exists in SRT_PRIVATE_CONTEXT
00292       std::string path = priv + "/Mad/data";
00293       void *dir_ptr = gSystem -> OpenDirectory(path.c_str());
00294 
00295       if(!dir_ptr){
00296         base=pub; // if it doesn't exist use SRT_PUBLIC
00297       }
00298       else base = priv;
00299     }
00300     else{
00301       base=pub;
00302     }
00303     if(base=="") {
00304       MSG("AnalysisInfoAna",Msg::kFatal)<<"No SRT_PUBLIC_CONTEXT set "
00305                                         <<"Do not know where to look "
00306                                         <<"for AB pdf files "<<std::endl;
00307       assert(false);
00308     }
00309     base+="/Mad/data";
00310     string fname=base;
00311     string rmc="";
00312     if(ReleaseType::IsCedar(release)&&ReleaseType::IsDaikon(release)){
00313       rmc="cedar_daikon";
00314     }
00315     else if(ReleaseType::IsCedar(release)&&ReleaseType::IsData(release)){
00316       rmc="cedar_daikon";
00317     }
00318     else{
00319       MSG("AnalysisInfoAna",Msg::kWarning)<<"Dont know reco/mc versions "
00320                                           <<"defaulting to cedar_daikon"<<endl;
00321       rmc="cedar_daikon";
00322     }
00323     string sdet="";
00324     if(fDetectorType==Detector::kNear){
00325       sdet="near";
00326     }
00327     else if(fDetectorType==Detector::kFar){
00328       sdet="far";
00329     }
00330     else{
00331       MSG("AnalysisInfoAna",Msg::kWarning)<<"Dont know detector type "
00332                                           <<"defaulting to far"<<endl;
00333       sdet="far";
00334     }
00335 
00336     string sbeam="";
00337     switch (beam){
00338     case BeamType::kL010z000i:     sbeam="le0";      break;
00339     case BeamType::kL010z170i:     sbeam="le170";    break;
00340     case BeamType::kL010z185i:     sbeam="le";       break;
00341     case BeamType::kL010z200i:     sbeam="le200";    break;
00342     case BeamType::kL100z200i:     sbeam="pme";      break;
00343     case BeamType::kL150z200i:     sbeam="pme";      break;
00344     case BeamType::kL250z200i:     sbeam="phe";      break;
00345     case BeamType::kLE:            sbeam="le";       break;
00346     default:
00347       MSG("AnalysisInfoAna",Msg::kWarning)<<"Don't know beam type "
00348                                           <<" defaulting to LE"<<endl;
00349       sbeam="le";
00350       break;
00351     }
00352                                                                                 
00353     fname+="/ab_pdf_"+sdet+"_"+sbeam+"_"+rmc+".root";
00354  
00355     return fname; 
00356 }
00357 
00358 string AnalysisInfoAna::BuildROPIDFile()
00359 {
00360     string base="";
00361     base=getenv("SRT_PRIVATE_CONTEXT");
00362     if(base!=""&&base!="."){
00363       // check if directory exists in SRT_PRIVATE_CONTEXT
00364       std::string path = base + "/Mad/data";
00365       void *dir_ptr = gSystem -> OpenDirectory(path.c_str());
00366       if(!dir_ptr){
00367         base=getenv("SRT_PUBLIC_CONTEXT"); // if not there use SRT_PUBLIC_CONTEXT
00368       }
00369     }
00370     else{
00371       base=getenv("SRT_PUBLIC_CONTEXT");
00372     }
00373     if(base=="") {
00374       MSG("AnalysisInfoAna",Msg::kFatal)<<"No SRT_PUBLIC_CONTEXT set "
00375                                         <<"Do not know where to look "
00376                                         <<"for RO pdf files "<<std::endl;
00377       assert(false);
00378     }
00379     base+="/Mad/data";
00380     //fname will depend on reco/mc version and detector, but not on beam
00381     string rmc="";
00382     if(ReleaseType::IsCedar(release)&&ReleaseType::IsDaikon(release)){
00383       rmc="cedar.daikon";
00384     }
00385     else if(ReleaseType::IsCedar(release)&&ReleaseType::IsData(release)){
00386       rmc="cedar.daikon";
00387     }
00388     else {
00389       MSG("AnalysisInfoAna",Msg::kWarning)<<"Dont know reco/mc versions "
00390                              <<"defaulting to cedar_daikon "<<endl;
00391       rmc="cedar.daikon";
00392     }
00393     string sdet="";
00394     if(fDetectorType==Detector::kNear){
00395       sdet="near";
00396     }
00397     else if(fDetectorType==Detector::kFar){
00398       sdet="far";
00399     }
00400     else{
00401       MSG("AnalysisInfoAna",Msg::kWarning)<<"Dont know detector type "
00402                                           <<"defaulting to far"<<endl;
00403       sdet="far";
00404     }
00405     string fname=base;
00406     fname+="/knn.train."+sdet+"."+rmc+".root";
00407                                                                                 
00408     return fname;
00409 }