loon: ANtpAnalysisInfoAna.cxx Source File

00001 
00017 #include "CandNtupleSR/NtpSRRecord.h"
00018 #include "CandNtupleSR/NtpSREvent.h"
00019 #include "MessageService/MsgService.h"
00020 #include "StandardNtuple/NtpStRecord.h"
00021 #include "NueAna/SntpHelpers.h"
00022 #include "NueAna/ANtpAnalysisInfoAna.h"
00023 #include "AnalysisNtuples/ANtpEventInfo.h"
00024 #include "AnalysisNtuples/ANtpTrackInfo.h"
00025 #include "AnalysisNtuples/ANtpShowerInfo.h"
00026 #include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
00027 #include "AnalysisNtuples/ANtpDefaultValue.h"
00028 #include "Mad/MadNsID.h"
00029 #include "Mad/MadDpID.h"
00030 
00031 MadNsID ANtpAnalysisInfoAna::nsid;
00032 MadDpID ANtpAnalysisInfoAna::dpid;
00033 
00034 
00035 CVSID("$ID:");
00036 
00037 ANtpAnalysisInfoAna::ANtpAnalysisInfoAna(ANtpAnalysisInfoNue &anai):
00038    fANtpAnalysisInfo(anai)
00039 {
00040     fDetectorType =  Detector::kUnknown;
00041 }
00042 
00043 ANtpAnalysisInfoAna::~ANtpAnalysisInfoAna()
00044 {}
00045 
00046 //void ANtpAnalysisInfoAna::Analyze(int evtn, NtpSRRecord *srobj, NtpMCRecord * /*mcobj*/, NtpTHRecord * /*thobj*/)
00047 void ANtpAnalysisInfoAna::Analyze(int evtn, RecRecordImp<RecCandHeader> *srobj)
00048 {
00049   MSG("ANtpAnalysisInfoAna",Msg::kDebug)<<"Entering ANtpAnalysisInfoAna::Analyze"<<endl;
00050           
00051         
00052   if(srobj==0){
00053       return;
00054   }  if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00055         ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00056      return;
00057   }
00058 //      Reset(srobj->GetHeader().GetSnarl(),evtn);
00059 
00060   NtpSREvent *event = 0;
00061   event = SntpHelpers::GetEvent(evtn,srobj);
00062   if(!event){
00063       MSG("ANtpAnalysisInfoAna",Msg::kError)<<"Couldn't get event "<<evtn
00064          <<" from Snarl "<<srobj->GetHeader().GetSnarl()<<endl;
00065      return;
00066   }
00067 
00068  const RecCandHeader *ntpHeader = &(srobj->GetHeader());
00069  fDetectorType = ntpHeader->GetVldContext().GetDetector();
00070 
00071  NtpSRTrack *track = 0;
00072  NtpSRShower *shower = 0; 
00073  NtpSREventSummary *eventSummary;
00074 
00075   //and now an ugly bit of code to deal with either NtpStRecord or NtpSRRecord
00076   bool foundst=false;
00077   NtpStRecord *st=dynamic_cast<NtpStRecord *>(srobj);
00078   ANtpInfoObjectFiller* filla =  new ANtpInfoObjectFiller();
00079   ANtpRecoNtpManipulator *ntpManipulator = 0;
00080 
00081   //ANtpRecoNtpManipulator ntpManipulator;
00082   NtpSRRecord *sr = 0; 
00083  
00084   if(st != 0)
00085   {
00086       foundst = true;
00087       ntpManipulator =  new ANtpRecoNtpManipulator(st);
00088       eventSummary = &(st->evthdr);
00089 //.SetRecord(st);
00090   }
00091   else{
00092       sr=dynamic_cast<NtpSRRecord *>(srobj);
00093       ntpManipulator =  new ANtpRecoNtpManipulator(sr);
00094       eventSummary = &(sr->evthdr);
00095 //.SetRecord(sr);
00096   }
00097  
00098   //instansiate a NtpHelper object to help you get the info you want
00099   filla->SetStripArray(ntpManipulator->GetStripArray());
00100     
00101   //set up which flags you want to use to determine the primary shower or track
00102   //a value of 0 for a flag means it will not be used
00103   ntpManipulator->SetPrimaryShowerCriteria(0,1); // nplanes, total pulse height
00104   ntpManipulator->SetPrimaryTrackCriteria(0,1,0); // nplanes, length, total pulse height           
00105   //get the primary shower for the event - if no track is present it
00106   //returns 0
00107   ANtpEventManipulator * ntpEventManipulator = 
00108                             ntpManipulator->GetEventManipulator();
00109  
00110   ntpEventManipulator->SetEventInSnarl(evtn);   
00111   event=ntpEventManipulator->GetEvent();
00112   shower = ntpEventManipulator->GetPrimaryShower();
00113   track = ntpEventManipulator->GetPrimaryTrack(); 
00114   
00115   FillNueAnalysisInformation(event, track, shower, &fANtpAnalysisInfo, filla);
00116 
00117 
00118   Double_t niki_cc_pid;
00119 
00120   
00121   if(fDetectorType==Detector::kFar ){
00122        dpid.SetPHCorrection(1.018);
00123   } 
00124 
00125   BeamType::BeamType_t current_beam = BeamType::kLE;
00126   if(dpid.ChoosePDFs(fDetectorType,current_beam))
00127     fANtpAnalysisInfo.dpCCPID = dpid.CalcPID(track, event, eventSummary, fDetectorType, 0);
00128 
00129  //  MadNsID nsid;
00130   if(nsid.ChooseWeightFile(fDetectorType,current_beam)){
00131        if(!nsid.GetPID(event,track,shower,st,fDetectorType,fANtpAnalysisInfo.nsCCPID))
00132            niki_cc_pid=-999;
00133    }
00134    else niki_cc_pid=-999;
00135 
00136 
00137   if(filla){
00138      delete filla;
00139      filla=0;
00140   }
00141   if(ntpManipulator){
00142      delete ntpManipulator;
00143      ntpManipulator = 0;
00144   }
00145 
00146   MSG("ANtpAnalysisInfoAna",Msg::kDebug)<<"Leaving ANtpAnalysisInfoAna::Analyze"<<endl;           
00147 }
00148 //----------------------------------------------------------------------
00149 void ANtpAnalysisInfoAna::FillNueAnalysisInformation(NtpSREvent *ntpEvent, NtpSRTrack *ntpTrack, NtpSRShower *ntpShower, ANtpAnalysisInfoNue *analysisInfoNue, ANtpInfoObjectFiller * filla)
00150 {
00151 //    ANtpInfoObjectFiller filla;
00152 //    ANtpHeaderInfo *antpHeader       = new ANtpHeaderInfo();
00153 //    ANtpTruthInfoBeam *antpTruth     = new ANtpTruthInfoBeam();
00154 //    ANtpAnalysisInfo *antpAnalysis   = new ANtpAnalysisInfo();
00155     ANtpEventInfo *antpEvent         = new ANtpEventInfo();
00156     ANtpTrackInfo *antpTrack         = new ANtpTrackInfo();
00157     ANtpShowerInfo *antpShower       = new ANtpShowerInfo();
00158 
00159     filla->FillEventInformation(ntpEvent, antpEvent);
00160     if(ntpTrack != 0)  filla->FillTrackInformation(ntpTrack,antpTrack);
00161     if(ntpShower != 0) filla->FillShowerInformation(ntpShower,antpShower);      
00162 
00163 
00164     analysisInfoNue->isFullyContained  = -10;
00165     analysisInfoNue->passesCuts        = 1;
00166     analysisInfoNue->recoEventLength   = TMath::Abs(antpEvent->endPlane -
00167                        antpEvent->begPlane);
00168 
00169     if(ntpTrack != 0){
00170      analysisInfoNue->recoTrackLength   = TMath::Abs(antpTrack->endPlane -
00171                        antpTrack->begPlane);
00172      analysisInfoNue->recoTrackMomentum = antpTrack->fitMomentum;
00173      analysisInfoNue->recoTrackRange    = antpTrack->rangeMomentum;
00174      analysisInfoNue->recoSigmaQoverP   = antpTrack->sigmaQoverP;
00175       
00176      analysisInfoNue->recoMuEnergy      = RecoMuEnergy(ntpTrack);
00177 
00178      if(ntpShower != 0){
00179        analysisInfoNue->recoShowerEnergy  = RecoShwEnergy(ntpShower);
00180        analysisInfoNue->recoNuEnergy      = ( analysisInfoNue->recoMuEnergy +
00181                         analysisInfoNue->recoShowerEnergy );
00182      }
00183      analysisInfoNue->recoQENuEnergy     = RecoQENuEnergy(ntpTrack);
00184      analysisInfoNue->recoQEQ2           = RecoQEQ2(ntpTrack);
00185 
00186      if (analysisInfoNue->recoNuEnergy>0) {
00187         analysisInfoNue->recoHadronicY    = ( analysisInfoNue->recoShowerEnergy /
00188                        analysisInfoNue->recoNuEnergy );
00189      }
00190         
00191        analysisInfoNue->recoMuDCosZVtx    = RecoMuDCosZ(ntpTrack);
00192        analysisInfoNue->recoNuDCos        = RecoMuDCosNeu(ntpTrack);
00193     }
00194 
00195     analysisInfoNue->recoVtxX          = antpEvent->vtxX;
00196     analysisInfoNue->recoVtxY          = antpEvent->vtxY;
00197     analysisInfoNue->recoVtxZ          = antpEvent->vtxZ;
00198 
00199     analysisInfoNue->isFullyContained  = 
00200          IsFidAll(antpEvent->vtxX, antpEvent->vtxY, antpEvent->vtxZ, ntpEvent);
00201 //       analysisInfoNue->isFullyContained  = 0;  
00202     //fiducial criteria on vtx
00203 
00204     analysisInfoNue->inFiducialVolume = IsFidVtxEvt(ntpEvent);
00205 //       analysisInfoNue->inFiducialVolume = 0;
00206     
00207 
00208     if(antpEvent){
00209       delete antpEvent;
00210       antpEvent=0;
00211     }
00212 
00213     if(antpTrack){
00214       delete antpTrack;
00215       antpTrack=0;
00216     }
00217     if(antpShower){
00218       delete antpShower;
00219       antpShower=0;
00220     }
00221   
00222     return;
00223 }
00224 
00225 //*******************************************************
00226 Float_t ANtpAnalysisInfoAna::RecoMuEnergy(NtpSRTrack *ntpTrack){
00227 
00228      if(IsFidAll(ntpTrack->vtx.x, ntpTrack->vtx.y, ntpTrack->vtx.z)
00229                     || ntpTrack->momentum.qp==0) {
00230          return sqrt(ntpTrack->momentum.range*ntpTrack->momentum.range
00231                   + 0.10555*0.10555);
00232      }
00233      else {
00234        if(ntpTrack->momentum.qp < 1e-5) return 10000.0;
00235        else 
00236           return sqrt(1./(ntpTrack->momentum.qp*ntpTrack->momentum.qp)
00237                                 + 0.10555*0.10555);
00238      }
00239 
00240      return 0;
00241 }
00242 
00243 //*******************************************************
00244 Float_t ANtpAnalysisInfoAna::RecoShwEnergy(NtpSRShower *ntpShower){
00245           //use SR reco
00246   Float_t theGeV=ntpShower->ph.gev;
00247   return theGeV;
00248 }
00249 
00250 //*******************************************************
00251 Float_t ANtpAnalysisInfoAna::RecoQENuEnergy(NtpSRTrack *ntpTrack){
00252         
00253         Float_t nucleonMass = 0.93956563; //mass of neutron by default
00254         if(GetChargeSign(ntpTrack)==1) nucleonMass = 0.93827231; //proton mass for nubar
00255         Float_t muonEnergy = RecoMuEnergy(ntpTrack);
00256         Float_t muonMass = 0.10555;
00257         if(muonEnergy < muonMass){
00258           MSG("ANtpAnalysisInfoAna",Msg::kError)
00259                  << "muon Energy < muon mass | Large reco failure"
00260                  << endl;
00261           return ANtpDefVal::kFloat;
00262         }
00263 
00264         if(TMath::Abs(muonEnergy) > 1e10){
00265            MSG("ANtpAnalysisInfoAna",Msg::kError)
00266                   << "muon Energy too big "
00267                   << muonEnergy<< " stopping"
00268                   << endl;
00269            return ANtpDefVal::kFloat;
00270         }
00271                 
00272         
00273         Float_t muonMomentum = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
00274         Float_t costhbl = RecoMuDCosNeu(ntpTrack);
00275 
00276         Float_t Eneu = ANtpDefVal::kFloat;
00277         if(nucleonMass - muonEnergy + muonMomentum*costhbl > 1e-8) 
00278         Eneu = (nucleonMass*muonEnergy - muonMass*muonMass/2.)
00279                 /(nucleonMass - muonEnergy + muonMomentum*costhbl);
00280         
00281         return Eneu;
00282 }
00283                                                          
00284 //********************************************************       
00285 Float_t ANtpAnalysisInfoAna::RecoQEQ2(NtpSRTrack *ntpTrack){
00286 
00287      Float_t Eneu = RecoQENuEnergy(ntpTrack);
00288      if(ANtpDefVal::IsDefault(Eneu)) return ANtpDefVal::kFloat;
00289      Float_t muonEnergy = RecoMuEnergy(ntpTrack);
00290      Float_t muonMass = 0.10555;
00291      Float_t muonMomentum = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
00292      Float_t costhbl = RecoMuDCosNeu(ntpTrack);
00293      return -2.*Eneu*(muonEnergy-muonMomentum*costhbl)+muonMass*muonMass;
00294   
00295 }
00296 
00297 //*******************************************************
00298 Int_t ANtpAnalysisInfoAna::GetChargeSign(NtpSRTrack *ntpTrack){
00299   if(RecoMuQP(ntpTrack)>0) return 1;
00300   return -1;
00301 }
00302 
00303 //*******************************************************
00304 Float_t ANtpAnalysisInfoAna::RecoMuQP(NtpSRTrack *ntpTrack){
00305   return ntpTrack->momentum.qp;
00306 }
00307 
00308 //*******************************************************
00309 Float_t ANtpAnalysisInfoAna::RecoMuDCosNeu(NtpSRTrack *ntpTrack){ //ds_mu/ds_neu
00310   Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam
00311   Float_t bl_y = sqrt(1. - bl_z*bl_z);
00312   Float_t costhbl = ntpTrack->vtx.dcosz*bl_z + ntpTrack->vtx.dcosy*bl_y;
00313         
00314   return  costhbl;
00315 }
00316 
00317 //*******************************************************
00318 Float_t ANtpAnalysisInfoAna::RecoMuDCosZ(NtpSRTrack *ntpTrack){ //dz/ds
00319    return ntpTrack->vtx.dcosz;
00320 }
00321 
00322 //*******************************************************
00323 Bool_t ANtpAnalysisInfoAna::IsFidVtx(NtpSRTrack *ntpTrack){
00324   if(fDetectorType == Detector::kFar){
00325           
00326      if(ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>29.4 ||   //ends
00327           (ntpTrack->vtx.z<16.5&&ntpTrack->vtx.z>14.5) ||  //between SMs
00328           sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00329                +(ntpTrack->vtx.y*ntpTrack->vtx.y))>3.5 ||
00330              sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00331                  +(ntpTrack->vtx.y*ntpTrack->vtx.y))<0.4) return false;
00332    
00333   }
00334   else if(fDetectorType == Detector::kNear){
00335 
00336       if(ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>6.5 ||
00337            sqrt(((ntpTrack->vtx.x-1.3)*(ntpTrack->vtx.x-1.3)) +
00338               (ntpTrack->vtx.y*ntpTrack->vtx.y))>1) return false;
00339        }                                                                                                                    
00340   return true;
00341 }
00342 
00343 //*******************************************************
00344 Int_t ANtpAnalysisInfoAna::IsFidVtxEvt(NtpSREvent *ntpEvent){
00345 
00346   Bool_t contained = true;
00347         
00348   Int_t test = 0;
00349   MSG("ANtpAnalysisInfoAna",Msg::kDebug) << "DetectorType " << fDetectorType << endl;
00350   if(fDetectorType == Detector::kNear)
00351        test = InsideNearFiducial(ntpEvent->vtx.x, ntpEvent->vtx.y,
00352                                   ntpEvent->vtx.z);
00353   if(fDetectorType == Detector::kFar)
00354        test = InsideFarFiducial(ntpEvent->vtx.x, ntpEvent->vtx.y,
00355                                   ntpEvent->vtx.z);
00356   if(test <= 0) contained = false;
00357   MSG("ANtpAnalysisInfoAna",Msg::kDebug) << " IsFidVtxEvt " << test << endl;
00358   return test;
00359 }
00360 
00361 Int_t ANtpAnalysisInfoAna::IsFidVtxEvt(NtpSREvent *ntpEvent, Int_t detType){
00362 
00363   Bool_t contained = true;
00364         
00365   Int_t test = 0;
00366   MSG("ANtpAnalysisInfoAna",Msg::kDebug) << "DetectorType " << detType << endl;
00367   if(detType == Detector::kNear)
00368        test = InsideNearFiducial(ntpEvent->vtx.x, ntpEvent->vtx.y,
00369                                   ntpEvent->vtx.z);
00370   if(detType == Detector::kFar)
00371        test = InsideFarFiducial(ntpEvent->vtx.x, ntpEvent->vtx.y,
00372                                   ntpEvent->vtx.z);
00373   if(test <= 0) contained = false;
00374   
00375   MSG("ANtpAnalysisInfoAna",Msg::kDebug) << " IsFidVtxEvt " << test << endl;
00376   return test;
00377 }
00378 
00379 void ANtpAnalysisInfoAna::Set3DHit(DeqFloat_t &x
00380                              , DeqFloat_t &y
00381                              , DeqFloat_t &z
00382                              , DeqFloat_t &e){
00383         
00384      if(x.size()!=0 && y.size()!=0 && z.size()!=0 && e.size()!=0){
00385  
00386            fX=x;
00387            fY=y;
00388            fZ=z;
00389            fE=e;
00390      }
00391 }
00392 
00393 
00394 Int_t ANtpAnalysisInfoAna::IsFidAll(Float_t vtxX, Float_t vtxY, Float_t vtxZ, NtpSREvent *event)
00395 {
00396   
00397   Bool_t contained = true;
00398   Int_t test = 0;
00399   if(fX.size()==0 || fY.size()==0 || fZ.size()==0 || fE.size()==0){
00400 //       MAXMSG("ANtpAnalysisInfoAna",Msg::kWarning, 1)
00401 //            << "3D Hits not set for event,"
00402 //            << "cannot evaluate isFullyContained."
00403 //            << endl;
00404         return ANtpDefVal::kInt;
00405     }           
00406 
00407   for(UInt_t i = 0; i < fX.size() && contained; i++){
00408      test = 0;
00409      Float_t x = fX[i] + vtxX;
00410      Float_t y = fY[i] + vtxY;
00411      Float_t z = fZ[i] + vtxZ;
00412      if(fDetectorType == Detector::kNear)
00413              test = InsideNearFiducial(x, y, z);
00414      if(fDetectorType == Detector::kFar)
00415              test = InsideFarFiducial(x, y, z);
00416      if(test <= 0) contained = false;
00417   }
00418 
00419    if(event) 
00420    if(contained && event->plane.n > 66) //this is approximately 4 meteres which is the 3D hit cutoff
00421    {
00422      if(fDetectorType == Detector::kNear)
00423                 test = InsideNearFiducial(event->end.x, event->end.y, event->end.z);
00424      if(fDetectorType == Detector::kFar)
00425                 test = InsideFarFiducial(event->end.x, event->end.y, event->end.z);
00426      if(test <= 0) contained = false;
00427    }      
00428   
00429   if(contained) test = 1;
00430   return test;
00431 }
00432      
00433 Int_t ANtpAnalysisInfoAna::InsideFarFiducial(Float_t x, Float_t y, Float_t z)
00434 {
00435   Float_t SuperModule1Beg =  0.35;
00436   Float_t SuperModule2Beg = 16.20;
00437   Float_t SuperModule1End = 14.57;
00438   Float_t SuperModule2End = 29.62;
00439 
00440   Float_t radialInner = 0.40;
00441   Float_t radialOuter = 3.87;
00442   Bool_t zContained = false;
00443   Bool_t xyContained = false;
00444 
00445   Float_t r = TMath::Sqrt(x*x + y*y);
00446     
00447   if( (z >= SuperModule1Beg && z <=SuperModule1End) || 
00448       (z >= SuperModule2Beg && z <=SuperModule2End) )
00449      zContained = true;
00450 
00451   if( r >= radialInner && r <= radialOuter)
00452      xyContained = true;
00453   
00454   Int_t retVal = 0;
00455   if(zContained && xyContained) retVal = 1;
00456   if(!zContained) retVal = -1;
00457   if(!xyContained) retVal -= 2;
00458  
00459   return retVal;  //  1 contained, -1 out of bounds z 
00460                   //  -2 oob xy, -3 oob both
00461 }
00462 
00463 Int_t ANtpAnalysisInfoAna::InsideNearFiducial(Float_t x, Float_t y, Float_t z)
00464 {
00465   Float_t SuperModule1Beg = 0.40;
00466   Float_t SuperModule1End = 6.50;
00467                                                                                 
00468   Float_t radialInner = 0;
00469   Float_t radialOuter = 1;
00470   Float_t xCenter = 1.4885;
00471   Float_t yCenter = 0.1397;
00472   Bool_t zContained = false;
00473   Bool_t xyContained = false;
00474 
00475   Float_t r = TMath::Sqrt((x-xCenter)*(x-xCenter) + (y-yCenter)*(y-yCenter));
00476   if( z >= SuperModule1Beg && z <=SuperModule1End)
00477      zContained = true;
00478                                                                                 
00479   if( r >= radialInner && r <= radialOuter)
00480      xyContained = true;
00481 
00482   Int_t retVal = 0;
00483   if(zContained && xyContained) retVal = 1;
00484   if(!zContained) retVal = -1;
00485   if(!xyContained) retVal -= 2;           
00486   
00487   return retVal;
00488 }
00489