WCSimVertexNtuple.cc

Go to the documentation of this file.
00001 #include "WCSimVertexNtuple.hh"
00002 #include "WCSimVertexGeometry.hh"
00003 #include "WCSimVertexFinder.hh"
00004 
00005 #include "TDirectory.h"
00006 #include "TFile.h"
00007 #include "TTree.h"
00008 #include "TMath.h"
00009 
00010 #include <iostream>
00011 #include <cmath>
00012 
00013 ClassImp(WCSimVertexNtuple)
00014 
00015 WCSimVertexNtuple::WCSimVertexNtuple() : WCSimNtuple(), 
00016   fRecoFile(0),
00017   fRecoTree(0)
00018 {
00019   
00020 }
00021   
00022 WCSimVertexNtuple::~WCSimVertexNtuple()
00023 {
00024  
00025 }
00026 
00027 void WCSimVertexNtuple::Reset()
00028 {
00029   std::cout << " *** WCSimVertexNtuple::Reset() *** " << std::endl;
00030 
00031   fRecoFile = 0;
00032   fRecoTree = 0;
00033 
00034   return;
00035 }
00036 
00037 void WCSimVertexNtuple::Fill( WCSimTrueEvent* trueEvent, WCSimRecoEvent* recoEvent )
00038 {
00039   this->ResetVariables();
00040   this->WriteVariables( trueEvent, recoEvent );
00041 }
00042   
00043 void WCSimVertexNtuple::ResetVariables()
00044 {
00045   // reset variables
00046   // ===============
00047   fRunNum = -1;
00048   fEventNum = -1;
00049   fTriggerNum = -1;
00050  
00051   fVtxX = 0.0;
00052   fVtxY = 0.0;
00053   fVtxZ = 0.0;
00054   fVtxTime = 0.0;
00055   fDirX = 0.0;
00056   fDirY = 0.0;
00057   fDirZ = 0.0;
00058 
00059   fTrueVtxX = 0.0;
00060   fTrueVtxY = 0.0;
00061   fTrueVtxZ = 0.0;
00062   fTrueDirX = 0.0;
00063   fTrueDirY = 0.0;
00064   fTrueDirZ = 0.0;
00065   fTrueLength = 0.0;
00066  
00067   index = -1;
00068   last = -1;
00069 
00070   fIsFiltered = 0;
00071 
00072   fDigitX = 0.0;
00073   fDigitY = 0.0;
00074   fDigitZ = 0.0;
00075   fDigitT = 0.0;
00076   fDigitQ = 0.0;
00077 
00078   fConeAngle = 0.0;
00079   fAngle = 0.0;
00080   fZenith = 0.0;
00081   fAzimuth = 0.0;
00082   fSolidAngle = 0.0;  
00083 
00084   fDistPoint = 0.0;
00085   fDistTrack = 0.0;
00086   fDistPhoton = 0.0;
00087   fDistScatter = 0.0;
00088 
00089   fDeltaTime = 0.0;
00090   fDeltaSigma = 0.0;
00091 
00092   fDeltaAngle = 0.0;
00093   fDeltaPoint = 0.0;
00094   fDeltaTrack = 0.0;
00095   fDeltaPhoton = 0.0;
00096   fDeltaScatter = 0.0;
00097 
00098   fPointPath = 0.0;
00099   fExtendedPath = 0.0;
00100   fPointResidual = 0.0;
00101   fExtendedResidual = 0.0;
00102 
00103   // --- debug ---
00104   fExtendedResidualCorrected = 0.0;
00105 
00106   return;
00107 }
00108 
00109 void WCSimVertexNtuple::WriteVariables( WCSimTrueEvent* fTrueEvent, WCSimRecoEvent* fRecoEvent )
00110 {
00111   std::cout << " *** WCSimVertexNtuple::WriteVariables() *** " << std::endl;
00112   
00113   // run vertex finder
00114   // =================
00115   // WCSimVertexFinder::Instance()->Reset();
00116   // WCSimRecoVertex* fRecoVertex = WCSimVertexFinder::Instance()->Run(fRecoEvent,
00117   //                                                                  fTrueEvent->GetVtxX(),
00118   //                                                                  fTrueEvent->GetVtxY(),
00119   //                                                                  fTrueEvent->GetVtxZ(),
00120   //                                                                  fTrueEvent->GetDirX(),
00121   //                                                                  fTrueEvent->GetDirY(),
00122   //                                                                  fTrueEvent->GetDirZ());
00123   //
00124 
00125   // get true information
00126   // ====================
00127   
00128   fTrueVtxX = fTrueEvent->GetVtxX();
00129   fTrueVtxY = fTrueEvent->GetVtxY();
00130   fTrueVtxZ = fTrueEvent->GetVtxZ();
00131   fTrueDirX = fTrueEvent->GetDirX();
00132   fTrueDirY = fTrueEvent->GetDirY();
00133   fTrueDirZ = fTrueEvent->GetDirZ();
00134   fTrueLength = fTrueEvent->GetLength();
00135 
00136   // get reconstructed vertex
00137   // ========================
00138   WCSimRecoVertex* fRecoVertex = (WCSimRecoVertex*)(fRecoEvent->GetVertex());
00139 
00140   // check vertex
00141   // ============
00142   if( fRecoVertex==0 ) return;
00143 
00144   // calculate vertex residuals
00145   // ==========================
00146   WCSimVertexGeometry::Instance()->CalcResiduals(fRecoEvent,fRecoVertex);
00147     
00148   // header information
00149   // ==================
00150   fRunNum = fRecoEvent->GetRun();
00151   fEventNum = fRecoEvent->GetEvent();
00152   fTriggerNum = fRecoEvent->GetTrigger();
00153 
00154   // get vertex position
00155   // ===================
00156   fVtxX = fRecoVertex->GetX();
00157   fVtxY = fRecoVertex->GetY();
00158   fVtxZ = fRecoVertex->GetZ();
00159   fVtxTime = fRecoVertex->GetTime();
00160   
00161   fDirX = fRecoVertex->GetDirX();
00162   fDirY = fRecoVertex->GetDirY();
00163   fDirZ = fRecoVertex->GetDirZ();
00164 
00165   // get digits used in fit
00166   // ======================
00167   last = WCSimVertexGeometry::Instance()->GetNDigits();
00168 
00169   for( Int_t idigit=0; idigit<last; idigit++ ){
00170     index = idigit;
00171 
00172     fIsFiltered = WCSimVertexGeometry::Instance()->IsFiltered(idigit);
00173 
00174     fDigitX = WCSimVertexGeometry::Instance()->GetDigitX(idigit);
00175     fDigitY = WCSimVertexGeometry::Instance()->GetDigitY(idigit);
00176     fDigitZ = WCSimVertexGeometry::Instance()->GetDigitZ(idigit);
00177     fDigitT = WCSimVertexGeometry::Instance()->GetDigitT(idigit);
00178     fDigitQ = WCSimVertexGeometry::Instance()->GetDigitQ(idigit);  
00179 
00180     fConeAngle   = WCSimVertexGeometry::Instance()->GetConeAngle(idigit);
00181     fAngle       = WCSimVertexGeometry::Instance()->GetAngle(idigit);
00182     fZenith      = WCSimVertexGeometry::Instance()->GetZenith(idigit);
00183     fAzimuth     = WCSimVertexGeometry::Instance()->GetAzimuth(idigit);
00184     fSolidAngle  = WCSimVertexGeometry::Instance()->GetSolidAngle(idigit);
00185 
00186     fDistPoint   = WCSimVertexGeometry::Instance()->GetDistPoint(idigit);
00187     fDistTrack   = WCSimVertexGeometry::Instance()->GetDistTrack(idigit);
00188     fDistPhoton  = WCSimVertexGeometry::Instance()->GetDistPhoton(idigit);
00189     fDistScatter = WCSimVertexGeometry::Instance()->GetDistScatter(idigit);
00190 
00191     fDeltaTime   = WCSimVertexGeometry::Instance()->GetDeltaTime(idigit);
00192     fDeltaSigma  = WCSimVertexGeometry::Instance()->GetDeltaSigma(idigit);
00193     
00194     fDeltaAngle  = WCSimVertexGeometry::Instance()->GetDeltaAngle(idigit);
00195     fDeltaPoint  = WCSimVertexGeometry::Instance()->GetDeltaPoint(idigit);
00196     fDeltaTrack  = WCSimVertexGeometry::Instance()->GetDeltaTrack(idigit);
00197     fDeltaPhoton = WCSimVertexGeometry::Instance()->GetDeltaPhoton(idigit);
00198     fDeltaScatter = WCSimVertexGeometry::Instance()->GetDeltaScatter(idigit);
00199     
00200     fPointPath = WCSimVertexGeometry::Instance()->GetPointPath(idigit);
00201     fExtendedPath = WCSimVertexGeometry::Instance()->GetExtendedPath(idigit);
00202     fPointResidual = WCSimVertexGeometry::Instance()->GetPointResidual(idigit);
00203     fExtendedResidual = WCSimVertexGeometry::Instance()->GetExtendedResidual(idigit);
00204 
00205     // --- debug ---
00206     fExtendedResidualCorrected = fExtendedResidual 
00207                                + WCSimVertexGeometry::Instance()->GetDeltaCorrection(idigit,fTrueLength);
00208 
00209     this->WriteToFile();
00210   }
00211 
00212   return;
00213 }
00214 
00215 void WCSimVertexNtuple::WriteToFile()
00216 {
00217   TDirectory* tmpd = 0;
00218 
00219   if( fRecoFile==0 ){
00220     tmpd = gDirectory;
00221     std::cout << " *** WCSimVertexNtuple::OpenFile(...) *** " << std::endl;
00222     std::cout << "  opening file: " << fNtpFileName.Data() << std::endl;
00223     fRecoFile = new TFile(fNtpFileName.Data(),"recreate");
00224     fRecoTree = new TTree("ntuple","my analysis ntuple");
00225     fRecoTree->Branch("run",&fRunNum,"run/I");
00226     fRecoTree->Branch("event",&fEventNum,"event/I");
00227     fRecoTree->Branch("trigger",&fTriggerNum,"trigger/I");
00228     fRecoTree->Branch("vtxX",&fVtxX,"vtxX/D");
00229     fRecoTree->Branch("vtxY",&fVtxY,"vtxY/D");
00230     fRecoTree->Branch("vtxZ",&fVtxZ,"vtxZ/D");
00231     fRecoTree->Branch("vtxTime",&fVtxTime,"vtxTime/D");
00232     fRecoTree->Branch("dirX",&fDirX,"dirX/D");
00233     fRecoTree->Branch("dirY",&fDirY,"dirY/D");
00234     fRecoTree->Branch("dirZ",&fDirZ,"dirZ/D");
00235     fRecoTree->Branch("trueVtxX",&fTrueVtxX,"trueVtxX/D");
00236     fRecoTree->Branch("trueVtxY",&fTrueVtxY,"trueVtxY/D");
00237     fRecoTree->Branch("trueVtxZ",&fTrueVtxZ,"trueVtxZ/D");
00238     fRecoTree->Branch("trueDirX",&fTrueDirX,"trueDirX/D");
00239     fRecoTree->Branch("trueDirY",&fTrueDirY,"trueDirY/D");
00240     fRecoTree->Branch("trueDirZ",&fTrueDirZ,"trueDirZ/D");
00241     fRecoTree->Branch("trueLength",&fTrueLength,"trueLength/D");
00242     fRecoTree->Branch("index",&index,"index/I");
00243     fRecoTree->Branch("last",&last,"last/I");
00244     fRecoTree->Branch("filter",&fIsFiltered,"filter/I");
00245     fRecoTree->Branch("digitX",&fDigitX,"digitX/D");
00246     fRecoTree->Branch("digitY",&fDigitY,"digitY/D");
00247     fRecoTree->Branch("digitZ",&fDigitZ,"digitZ/D");
00248     fRecoTree->Branch("digitT",&fDigitT,"digitT/D");
00249     fRecoTree->Branch("digitQ",&fDigitQ,"digitQ/D");
00250     fRecoTree->Branch("coneAngle",&fConeAngle,"coneAngle/D");
00251     fRecoTree->Branch("angle",&fAngle,"angle/D");
00252     fRecoTree->Branch("zenith",&fZenith,"zenith/D");
00253     fRecoTree->Branch("azimuth",&fAzimuth,"azimuth/D");
00254     fRecoTree->Branch("solidAngle",&fSolidAngle,"solidAngle/D");
00255     fRecoTree->Branch("distPoint",&fDistPoint,"distPoint/D");
00256     fRecoTree->Branch("distTrack",&fDistTrack,"distTrack/D");
00257     fRecoTree->Branch("distPhoton",&fDistPhoton,"distPhoton/D");
00258     fRecoTree->Branch("distScatter",&fDistScatter,"distScatter/D");
00259     fRecoTree->Branch("deltaTime",&fDeltaTime,"deltaTime/D");
00260     fRecoTree->Branch("deltaSigma",&fDeltaSigma,"deltaSigma/D");
00261     fRecoTree->Branch("deltaAngle",&fDeltaAngle,"deltaAngle/D");
00262     fRecoTree->Branch("deltaPoint",&fDeltaPoint,"deltaPoint/D");
00263     fRecoTree->Branch("deltaTrack",&fDeltaTrack,"deltaTrack/D");
00264     fRecoTree->Branch("deltaPhoton",&fDeltaPhoton,"deltaPhoton/D");
00265     fRecoTree->Branch("deltaScatter",&fDeltaScatter,"deltaScatter/D");
00266     fRecoTree->Branch("pointPath",&fPointPath,"pointPath/D");
00267     fRecoTree->Branch("extendedPath",&fExtendedPath,"extendedPath/D"); 
00268     fRecoTree->Branch("pointResidual",&fPointResidual,"pointResidual/D");
00269     fRecoTree->Branch("extendedResidual",&fExtendedResidual,"extendedResidual/D"); 
00270 
00271     // --- debug ---
00272     fRecoTree->Branch("extendedResidualCorrected",&fExtendedResidualCorrected,"extendedResidualCorrected/D"); 
00273 
00274     gDirectory = tmpd;
00275   }
00276 
00277   if( fRecoFile ){
00278     tmpd = gDirectory;
00279     fRecoFile->cd();
00280     fRecoTree->Fill();
00281     gDirectory = tmpd;
00282   }
00283 
00284   return;
00285 }
00286 
00287 void WCSimVertexNtuple::CloseFile()
00288 {
00289   TDirectory* tmpd = 0;
00290 
00291   if( fRecoFile ){
00292     tmpd = gDirectory;
00293     std::cout << " *** WCSimVertexNtuple::CloseFile(...) *** " << std::endl;
00294     std::cout << "  closing file: " << fNtpFileName.Data() << std::endl;
00295     fRecoFile->cd();
00296     fRecoTree->Write();
00297     fRecoFile->Close();
00298     gDirectory = tmpd;
00299   }
00300 
00301   return;
00302 }