loon: VtxFinderAna.cxx Source File

00001 #include "TH1F.h"
00002 #include "TF1.h"
00003 #include "TCanvas.h"
00004 #include "TClonesArray.h"
00005 #include "TMinuit.h"
00006 #include "StandardNtuple/NtpStRecord.h"
00007 #include "CandNtupleSR/NtpSRRecord.h"
00008 #include "CandNtupleSR/NtpSREvent.h"
00009 #include "CandNtupleSR/NtpSRTrack.h"
00010 #include "CandNtupleSR/NtpSRStrip.h"
00011 #include "MessageService/MsgService.h"
00012 #include "NueAna/NueAnaTools/SntpHelpers.h"
00013 #include "AnalysisNtuples/ANtpDefaultValue.h"
00014 #include "MCNtuple/NtpMCTruth.h"
00015 #include "MCNtuple/NtpMCStdHep.h"
00016 #include "TruthHelperNtuple/NtpTHEvent.h"
00017 #include "TruthHelperNtuple/NtpTHTrack.h"
00018 #include "TruthHelperNtuple/NtpTHShower.h"
00019 #include "TruthHelperNtuple/NtpTHStrip.h"
00020 #include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
00021 
00022 #include "NueAna/AnalysisInfoAna.h"
00023 #include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
00024 #include "VertexFinder/Module/VtxFinderAna.h"
00025                         
00026 CVSID("$Id: VtxFinderAna.cxx,v 1.4 2008/07/11 19:19:48 boehm Exp $");
00027 
00028 //  Retrieve teh event information
00029 //   get the right trk if any 
00030 //   get the mc information (if any)                             
00031 
00032 VtxFinderAna::VtxFinderAna(VtxFinderData &vf):
00033   fVtxData(vf)
00034 {
00035 }
00036 
00037 VtxFinderAna::~VtxFinderAna()
00038 {
00039 }
00040 
00041 void VtxFinderAna::Analyze(int event, RecRecordImp<RecCandHeader> *srobj)
00042 {
00043   if(srobj==0){
00044     return;
00045   }
00046 
00047   NtpStRecord *st = dynamic_cast<NtpStRecord *>(srobj);
00048   vtxf = new NtpVtxFinder(event, st);
00049 
00050   Analyze(event,st, srobj);
00051   delete vtxf;
00052 }
00053 
00054 
00055 void VtxFinderAna::Analyze(int evtn, NtpStRecord *srobj, RecRecordImp<RecCandHeader> *srobj2)
00056 {
00057   if(srobj==0){
00058     return;
00059   }
00060 
00061    NtpMCTruth *mctruth = 0;
00062    NtpMCStdHep *mcstdhep = 0;
00063    NtpTHEvent *thevent = 0;
00064    ANtpRecoNtpManipulator ntpManipulator(srobj);
00065 
00066    thevent = ntpManipulator.GetMCManipulator()->GetNtpTHEvent(evtn);
00067    if(thevent){
00068       mctruth = ntpManipulator.GetMCManipulator()->GetNtpMCTruth(thevent->neumc);
00069       mcstdhep = ntpManipulator.GetMCManipulator()->GetNtpMCStdHep(thevent->neustdhep);
00070     }
00071 
00072     NtpSREvent *event = 0;
00073     NtpSRTrack *track = 0;
00074                                                                                           
00075     //set up which flags you want to use to determine the primary shower or track
00076     //a value of 0 for a flag means it will not be used
00077     ntpManipulator.SetPrimaryTrackCriteria(0,1,0); // nplanes, length, total pulse height
00078     //get the primary track for the event - if no track is present it
00079                        //                        //returns 0
00080     ANtpEventManipulator * ntpEventManipulator =
00081                 ntpManipulator.GetEventManipulator();
00082                                                                                           
00083      ntpEventManipulator->SetEventInSnarl(evtn);
00084      event = ntpEventManipulator->GetEvent();
00085      track = ntpEventManipulator->GetPrimaryTrack();
00086 
00087 //In theory all relevant quantities have now been loaded
00088               
00089     if(vtxf->FindVertex() < 0) return;
00090      
00091     Float_t x = 0;
00092     Float_t y = 0; 
00093                                   
00094     fVtxData.vfVtxPlane = vtxf->VtxPlane();
00095     fVtxData.vfVtxX = vtxf->VtxX();                     
00096     fVtxData.vfVtxY = vtxf->VtxY();
00097     fVtxData.vfVtxZ = vtxf->VtxZ();
00098     fVtxData.vfVtxU = vtxf->VtxU();
00099     fVtxData.vfVtxV = vtxf->VtxV();
00100 
00101     if(track){
00102      fVtxData.trkVtxX = x = track->vtx.x;
00103      fVtxData.trkVtxY = y = track->vtx.y;
00104      fVtxData.trkVtxZ = track->vtx.z;
00105      fVtxData.trkVtxU = (sqrt(2.0)/2.0)*(x + y);
00106      fVtxData.trkVtxV = (sqrt(2.0)/2.0)*(-x + y);
00107     }
00108 
00109     if(mctruth){
00110       fVtxData.mcVtxX = x = mctruth->vtxx;
00111       fVtxData.mcVtxY = y = mctruth->vtxy;
00112       fVtxData.mcVtxZ = mctruth->vtxz;
00113       fVtxData.mcVtxU = (sqrt(2.0)/2.0)*(x + y);
00114       fVtxData.mcVtxV = (sqrt(2.0)/2.0)*(-x + y); 
00115 
00116       Float_t dZ = fVtxData.vfVtxZ - fVtxData.mcVtxZ;
00117       Float_t dU = fVtxData.vfVtxU - fVtxData.mcVtxU;
00118       Float_t dV = fVtxData.vfVtxV - fVtxData.mcVtxV;
00119 
00120       fVtxData.vfErrZ = dZ;
00121       fVtxData.vfErrU = dU;
00122       fVtxData.vfErrV = dV;                                                     
00123       fVtxData.vfError = TMath::Sqrt(dZ*dZ + dU*dU + dV*dV);
00124 
00125       FindClosest(event, srobj2);
00126       
00127       if(track){
00128         dZ = fVtxData.trkVtxZ - fVtxData.mcVtxZ;
00129         Float_t dY = fVtxData.trkVtxY - fVtxData.mcVtxY;
00130         Float_t dX = fVtxData.trkVtxX - fVtxData.mcVtxX;
00131                                                                                 
00132         fVtxData.trkError = TMath::Sqrt(dZ*dZ + dY*dY + dX*dX);
00133       }
00134     }
00135 
00136                                                                                 
00137     const RecCandHeader *ntpHeader = &(srobj2->GetHeader());
00138     Int_t detType = ntpHeader->GetVldContext().GetDetector();
00139 
00140     if(detType == Detector::kNear){
00141        fVtxData.mcContained = InsideNearFiducial(fVtxData.mcVtxX, fVtxData.mcVtxY, fVtxData.mcVtxZ);
00142        fVtxData.vfContained = InsideNearFiducial(fVtxData.vfVtxX, fVtxData.vfVtxY, fVtxData.vfVtxZ);
00143     }
00144     if(detType == Detector::kFar){
00145        fVtxData.mcContained = InsideFarFiducial(fVtxData.mcVtxX, fVtxData.mcVtxY, fVtxData.mcVtxZ);
00146        fVtxData.vfContained = InsideFarFiducial(fVtxData.vfVtxX, fVtxData.vfVtxY, fVtxData.vfVtxZ);
00147     }
00148 
00149 }
00150 
00151 
00152 void VtxFinderAna::FindClosest(NtpSREvent* event, 
00153                                   RecRecordImp<RecCandHeader> *srobj)
00154 {
00155   Double_t dist = 100000.0;
00156   fVtxData.vfBestU = ANtpDefaultValue::kFloat;
00157   fVtxData.vfBestV = ANtpDefaultValue::kFloat;
00158   fVtxData.vfBestZ = ANtpDefaultValue::kFloat;
00159                                                                                                          
00160   for(int i=0;i<event->nstrip;i++){
00161       Int_t index = SntpHelpers::GetStripIndex(i,event);
00162       NtpSRStrip *strip = SntpHelpers::GetStrip(index,srobj);
00163                                                                                       Double_t zpos = strip->z;
00164       Float_t dZ = TMath::Abs(zpos - fVtxData.mcVtxZ);
00165       if(dZ < dist){
00166          fVtxData.vfBestZ = zpos;
00167          dist = dZ;
00168       }                                                                        
00169   }
00170 
00171   Float_t utemp = 100000.0;
00172   Float_t vtemp = utemp;
00173   for(int i=0;i<event->nstrip;i++){
00174       Int_t index = SntpHelpers::GetStripIndex(i,event);
00175       NtpSRStrip *strip = SntpHelpers::GetStrip(index,srobj);
00176 
00177      
00178      if(TMath::Abs(fVtxData.vfBestZ - strip->z) < 0.25){   
00179        Double_t tpos = strip->tpos;
00180 
00181        if(strip->planeview==PlaneView::kU){
00182           Float_t dU = TMath::Abs(tpos - fVtxData.mcVtxU);
00183           if(dU < utemp){
00184             fVtxData.vfBestU = tpos;
00185             utemp = dU;
00186           }
00187        } else if(strip->planeview==PlaneView::kV){
00188           Float_t dV = TMath::Abs(tpos - fVtxData.mcVtxV);
00189           MSG("VtxFinderAna", Msg::kVerbose)<<"V problem: "
00190                    <<"dV = "<<dV<<" tpos =  "<<tpos
00191                    <<" vtemp = "<<vtemp<<" mc.vtxY = "
00192                    <<fVtxData.mcVtxV<<endl;
00193           if(dV < vtemp){
00194             fVtxData.vfBestV = tpos;
00195             vtemp = dV;
00196           }
00197       } 
00198     }
00199   }
00200 
00201   Float_t dU = TMath::Abs(fVtxData.vfBestU - fVtxData.mcVtxU);
00202   Float_t dV = TMath::Abs(fVtxData.vfBestV - fVtxData.mcVtxV);
00203   Float_t dZ = TMath::Abs(fVtxData.vfBestZ - fVtxData.mcVtxZ);
00204 
00205   fVtxData.vfBestErrU = dU;
00206   fVtxData.vfBestErrV = dV;
00207   fVtxData.vfBestErrZ = dZ;
00208 
00209   dist = TMath::Sqrt(dZ*dZ + dU*dU + dV*dV);
00210   fVtxData.vfBestError = dist;
00211 
00212 
00213   dU = TMath::Abs(fVtxData.vfBestU - fVtxData.vfVtxU);
00214   dV = TMath::Abs(fVtxData.vfBestV - fVtxData.vfVtxV);
00215   dZ = TMath::Abs(fVtxData.vfBestZ - fVtxData.vfVtxZ);
00216 
00217   fVtxData.vfRelErrU = dU;
00218   fVtxData.vfRelErrV = dV;
00219   fVtxData.vfRelErrZ = dZ;
00220   fVtxData.vfRelError =  TMath::Sqrt(dZ*dZ + dU*dU + dV*dV);
00221 
00222 
00223   Float_t u = fVtxData.vfBestU;
00224   Float_t v = fVtxData.vfBestV;
00225 
00226   if(ANtpDefVal::IsDefault(fVtxData.vfBestU) || 
00227      ANtpDefVal::IsDefault(fVtxData.vfBestV) ||
00228         ANtpDefVal::IsDefault(fVtxData.vfBestZ) )
00229   {
00230      fVtxData.vfBestU = ANtpDefaultValue::kFloat;
00231      fVtxData.vfBestV = ANtpDefaultValue::kFloat;
00232      fVtxData.vfBestZ = ANtpDefaultValue::kFloat;
00233      fVtxData.vfBestError = ANtpDefaultValue::kFloat;
00234      return;
00235   }
00236 
00237   fVtxData.vfBestX = (sqrt(2.0)/2.0)*(u-v);
00238   fVtxData.vfBestY = (sqrt(2.0)/2.0)*(u+v);                                         
00239 }
00240 
00241 Int_t VtxFinderAna::InsideFarFiducial(Float_t x, Float_t y, Float_t z)
00242 {
00243   Float_t SuperModule1Beg =  0.35;
00244   Float_t SuperModule2Beg = 16.20;
00245   Float_t SuperModule1End = 14.57;
00246   Float_t SuperModule2End = 29.62;
00247                                                                                 
00248   Float_t radialInner = 0.40;
00249   Float_t radialOuter = 3.87;
00250   Bool_t zContained = false;
00251   Bool_t xyContained = false;
00252                                                                                 
00253   Float_t r = TMath::Sqrt(x*x + y*y);
00254                                                                                 
00255   if( (z >= SuperModule1Beg && z <=SuperModule1End) ||
00256       (z >= SuperModule2Beg && z <=SuperModule2End) )
00257      zContained = true;
00258                                                                                 
00259   if( r >= radialInner && r <= radialOuter)
00260      xyContained = true;
00261                                                                                 
00262   Int_t retVal = 0;
00263   if(zContained && xyContained) retVal = 1;
00264   if(!zContained) retVal = -1;
00265   if(!xyContained) retVal -= 2;
00266                                                                                 
00267   return retVal;  //  1 contained, -1 out of bounds z
00268                   //  -2 oob xy, -3 oob both
00269 }
00270                                                                                 
00271 Int_t VtxFinderAna::InsideNearFiducial(Float_t x, Float_t y, Float_t z)
00272 {
00273   Float_t SuperModule1Beg = 0.40;
00274   Float_t SuperModule1End = 6.50;
00275                                                                                 
00276   Float_t radialInner = 0;
00277   Float_t radialOuter = 1;
00278   Float_t xCenter = 1.4885;
00279   Float_t yCenter = 0.1397;
00280   Bool_t zContained = false;
00281   Bool_t xyContained = false;
00282                                                                                 
00283   Float_t r = TMath::Sqrt((x-xCenter)*(x-xCenter) + (y-yCenter)*(y-yCenter));
00284   if( z >= SuperModule1Beg && z <=SuperModule1End)
00285      zContained = true;
00286                                                                                 
00287   if( r >= radialInner && r <= radialOuter)
00288      xyContained = true;
00289                                                                                 
00290   Int_t retVal = 0;
00291   if(zContained && xyContained) retVal = 1;
00292   if(!zContained) retVal = -1;
00293   if(!xyContained) retVal -= 2;
00294                                                                                 
00295   return retVal;
00296 }
00297