WCSimVertexSeedNtuple.cc

Go to the documentation of this file.
00001 #include "WCSimVertexSeedNtuple.hh"
00002 #include "WCSimVertexGeometry.hh"
00003 
00004 #include "TDirectory.h"
00005 #include "TFile.h"
00006 #include "TTree.h"
00007 #include "TMath.h"
00008 
00009 #include <iostream>
00010 #include <cmath>
00011 
00012 ClassImp(WCSimVertexSeedNtuple)
00013 
00014 WCSimVertexSeedNtuple::WCSimVertexSeedNtuple() : WCSimNtuple(),
00015   fRecoFile(0),
00016   fRecoTree(0)
00017 {
00018 
00019 }
00020 
00021 WCSimVertexSeedNtuple::~WCSimVertexSeedNtuple()
00022 {
00023 
00024 }
00025 
00026 void WCSimVertexSeedNtuple::Reset()
00027 {
00028   std::cout << " *** WCSimVertexSeedNtuple::Reset() *** " << std::endl;
00029 
00030   fRecoFile = 0;
00031   fRecoTree = 0;
00032 
00033   return;
00034 }
00035 
00036 void WCSimVertexSeedNtuple::Fill( WCSimTrueEvent* trueEvent, WCSimRecoEvent* recoEvent )
00037 {
00038   this->ResetVariables();
00039   this->WriteVariables( trueEvent, recoEvent );
00040 }
00041 
00042 void WCSimVertexSeedNtuple::ResetVariables()
00043 {
00044  // reset variables
00045   // ===============
00046   fRunNum = -1;
00047   fEventNum = -1;
00048   fTriggerNum = -1;
00049 
00050   index = -1;
00051   last = -1;
00052 
00053   fTrueVtxX = 0.0;
00054   fTrueVtxY = 0.0;
00055   fTrueVtxZ = 0.0;
00056   fTrueDirX = 0.0;
00057   fTrueDirY = 0.0;
00058   fTrueDirZ = 0.0;
00059 
00060   fRecoVtxX = 0.0;
00061   fRecoVtxY = 0.0;
00062   fRecoVtxZ = 0.0;
00063   fRecoVtxTime = 0.0;
00064 
00065   fRecoDelta = -99999.9;
00066   fRecoDeltaTrans = -99999.9;
00067   fRecoDeltaLong = -99999.9;
00068 
00069   return;
00070 }
00071 
00072 void WCSimVertexSeedNtuple::WriteVariables( WCSimTrueEvent* fTrueEvent, WCSimRecoEvent* fRecoEvent )
00073 {
00074   std::cout << " *** WCSimVertexSeedNtuple::WriteVariables() *** " << std::endl;
00075 
00076   // number of seeds (try 1000)
00077   // ==========================
00078   Int_t fNumSeeds = 1000;
00079 
00080   // calculate vertex seeds
00081   // ======================
00082   WCSimVertexGeometry::Instance()->CalcVertexSeeds(fRecoEvent,fNumSeeds);
00083 
00084   // header information
00085   // ==================
00086   fRunNum = fRecoEvent->GetRun();
00087   fEventNum = fRecoEvent->GetEvent();
00088   fTriggerNum = fRecoEvent->GetTrigger();
00089 
00090   // true vertex position
00091   // ====================
00092   fTrueVtxX  = fTrueEvent->GetVtxX();
00093   fTrueVtxY  = fTrueEvent->GetVtxY();
00094   fTrueVtxZ  = fTrueEvent->GetVtxZ();
00095   fTrueDirX  = fTrueEvent->GetDirX();
00096   fTrueDirY  = fTrueEvent->GetDirY();
00097   fTrueDirZ  = fTrueEvent->GetDirZ();
00098 
00099   // loop over vertex seeds
00100   // ======================
00101   UInt_t nlast = WCSimVertexGeometry::Instance()->GetNSeeds();
00102 
00103   for( UInt_t n=0; n<nlast; n++ ){
00104     index = (Int_t)n;
00105     last = (Int_t)nlast;
00106 
00107     fRecoVtxX = WCSimVertexGeometry::Instance()->GetSeedVtxX(n);
00108     fRecoVtxY = WCSimVertexGeometry::Instance()->GetSeedVtxY(n);
00109     fRecoVtxZ = WCSimVertexGeometry::Instance()->GetSeedVtxZ(n);
00110     fRecoVtxTime = WCSimVertexGeometry::Instance()->GetSeedVtxTime(n);
00111 
00112     Double_t dx = fRecoVtxX - fTrueVtxX;
00113     Double_t dy = fRecoVtxY - fTrueVtxY;
00114     Double_t dz = fRecoVtxZ - fTrueVtxZ;
00115     Double_t ds = sqrt( dx*dx + dy*dy + dz*dz );
00116     Double_t px = fTrueDirX;
00117     Double_t py = fTrueDirY;
00118     Double_t pz = fTrueDirZ;
00119 
00120     Double_t epsilon = 1.0e-7;
00121 
00122     if( ds>epsilon ){
00123       px = dx/ds;
00124       py = dy/ds;
00125       pz = dz/ds;
00126     }
00127 
00128     Double_t costh = px*fTrueDirX+py*fTrueDirY+pz*fTrueDirZ;
00129     Double_t sinth = 0.0;
00130     if( costh<1.0-epsilon ) sinth = sqrt(1.0-costh*costh);
00131     else costh = 1.0;
00132 
00133     fRecoDelta = ds;
00134     fRecoDeltaTrans = ds*sinth;
00135     fRecoDeltaLong = ds*costh;
00136 
00137     this->WriteToFile();
00138   }
00139 
00140   return;
00141 }
00142 
00143 void WCSimVertexSeedNtuple::WriteToFile()
00144 {
00145   TDirectory* tmpd = 0;
00146 
00147   if( fRecoFile==0 ){
00148     tmpd = gDirectory;
00149     std::cout << " *** WCSimVertexSeedNtuple::OpenFile(...) *** " << std::endl;
00150     std::cout << "  opening file: " << fNtpFileName.Data() << std::endl;
00151     fRecoFile = new TFile(fNtpFileName.Data(),"recreate");
00152     fRecoTree = new TTree("ntuple","my analysis ntuple");
00153     fRecoTree->Branch("run",&fRunNum,"run/I");
00154     fRecoTree->Branch("event",&fEventNum,"event/I");
00155     fRecoTree->Branch("trigger",&fTriggerNum,"trigger/I");    
00156     fRecoTree->Branch("index",&index,"index/I");
00157     fRecoTree->Branch("last",&last,"last/I");
00158     fRecoTree->Branch("trueVtxX",&fTrueVtxX,"trueVtxX/D");
00159     fRecoTree->Branch("trueVtxY",&fTrueVtxY,"trueVtxY/D");
00160     fRecoTree->Branch("trueVtxZ",&fTrueVtxZ,"trueVtxZ/D");   
00161     fRecoTree->Branch("trueDirX",&fTrueDirX,"trueDirX/D");
00162     fRecoTree->Branch("trueDirY",&fTrueDirY,"trueDirY/D");
00163     fRecoTree->Branch("trueDirZ",&fTrueDirZ,"trueDirZ/D");
00164     fRecoTree->Branch("recoVtxX",&fRecoVtxX,"recoVtxX/D");
00165     fRecoTree->Branch("recoVtxY",&fRecoVtxY,"recoVtxY/D");
00166     fRecoTree->Branch("recoVtxZ",&fRecoVtxZ,"recoVtxZ/D"); 
00167     fRecoTree->Branch("recoVtxTime",&fRecoVtxTime,"recoVtxTime/D");  
00168     fRecoTree->Branch("recoDelta",&fRecoDelta,"recoDelta/D");
00169     fRecoTree->Branch("recoDeltaTrans",&fRecoDeltaTrans,"recoDeltaTrans/D");
00170     fRecoTree->Branch("recoDeltaLong",&fRecoDeltaLong,"recoDeltaLong/D");
00171     gDirectory = tmpd;
00172   }
00173 
00174   if( fRecoFile ){
00175     tmpd = gDirectory;
00176     fRecoFile->cd();
00177     fRecoTree->Fill();
00178     gDirectory = tmpd;
00179   }
00180 
00181   return;
00182 }
00183 
00184 void WCSimVertexSeedNtuple::CloseFile()
00185 {
00186   TDirectory* tmpd = 0;
00187 
00188   if( fRecoFile ){
00189     tmpd = gDirectory;
00190     std::cout << " *** WCSimVertexSeedNtuple::CloseFile() *** " << std::endl;
00191     std::cout << "  closing file: " << fNtpFileName.Data() << std::endl;
00192     fRecoFile->cd();
00193     fRecoTree->Write();
00194     fRecoFile->Close();
00195     gDirectory = tmpd;
00196   }
00197 
00198   return;
00199 }
00200