#include "ITruthVerticesModule.hxx" #include // Root #include // oaEvent #include #include #include // oaGeomInfo #include #include ClassImp(COMET::ITruthVerticesModule); ClassImp(COMET::ITruthVerticesModule::ITruthVertex); #define CVSTAG "\ $Name: v5r19 $" #define CVSID "\ $Id: ITruthVerticesModule.cxx,v 1.24 2012/01/25 07:00:33 marks Exp $" COMET::ITruthVerticesModule::ITruthVerticesModule(const char* name, const char* title) : fMaxNVertices(500), fNVtx(0), fNVtxFGD1(0), fNVtxFGD2(0), fNVtxP0D(0), fNVtxDsECal(0), fNVtxBrlECal(0), fNVtxP0DECal(0), fNVtxTPC1(0), fNVtxTPC2(0), fNVtxTPC3(0), fNVtxSMRD(0), fNVtxRestOfOffAxis(0), fNVtxINGRID(0), fNVtxOther(0), fVertices( new TClonesArray("COMET::ITruthVerticesModule::ITruthVertex", fMaxNVertices) ) { SetNameTitle(name, title); fIsEnabled = kTRUE; fDescription = "Module storing the truth information for primary vertices"; fCVSTagName = CVSTAG; fCVSID = CVSID; } COMET::ITruthVerticesModule::~ITruthVerticesModule() { delete fVertices; } void COMET::ITruthVerticesModule::InitializeBranches() { fOutputTree->Branch("NVtx", &fNVtx, "NVtx/I", fBufferSize); fOutputTree->Branch("NVtxFGD1", &fNVtxFGD1, "NVtxFGD1/I", fBufferSize); fOutputTree->Branch("NVtxFGD2", &fNVtxFGD2, "NVtxFGD2/I", fBufferSize); fOutputTree->Branch("NVtxP0D", &fNVtxP0D, "NVtxP0D/I", fBufferSize); fOutputTree->Branch("NVtxDsECal", &fNVtxDsECal, "NVtxDsECal/I", fBufferSize); fOutputTree->Branch("NVtxBrlECal", &fNVtxBrlECal, "NVtxBrlECal/I", fBufferSize); fOutputTree->Branch("NVtxP0DECal", &fNVtxP0DECal, "NVtxP0DECal/I", fBufferSize); fOutputTree->Branch("NVtxTPC1", &fNVtxTPC1, "NVtxTPC1/I", fBufferSize); fOutputTree->Branch("NVtxTPC2", &fNVtxTPC2, "NVtxTPC2/I", fBufferSize); fOutputTree->Branch("NVtxTPC3", &fNVtxTPC3, "NVtxTPC3/I", fBufferSize); fOutputTree->Branch("NVtxSMRD", &fNVtxSMRD, "NVtxSMRD/I", fBufferSize); fOutputTree->Branch("NVtxRestOfOffAxis", &fNVtxRestOfOffAxis, "NVtxRestOfOffAxis/I", fBufferSize); fOutputTree->Branch("NVtxINGRID", &fNVtxINGRID, "NVtxINGRID/I", fBufferSize); fOutputTree->Branch("NVtxOther", &fNVtxOther, "NVtxOther/I", fBufferSize); fOutputTree->Branch("Vertices", "TClonesArray", &fVertices, fBufferSize, fSplitLevel); } Bool_t COMET::ITruthVerticesModule::ProcessFirstEvent(COMET::ICOMETEvent& event) { bool IsDataFile = event.GetContext().IsDetector(); if(IsDataFile) { throw COMET::EDataFile(); } return true; } bool COMET::ITruthVerticesModule::FillTree(COMET::ICOMETEvent& event) { if( ! event.GetContext().IsMC() ) { COMETInfo("Event Context reports event is non-MC. Assuming entire file is non-MC."); return false; // disable module } if( ! event.Has("truth/G4PrimVertex00") ) { COMETWarn("No vertices container found, skipping event."); return true; // this happens occasionally, so don't disable the module } fNVtx = 0; fNVtxFGD1 = 0; fNVtxFGD2 = 0; fNVtxP0D = 0; fNVtxDsECal = 0; fNVtxBrlECal = 0; fNVtxP0DECal = 0; fNVtxTPC1 = 0; fNVtxTPC2 = 0; fNVtxTPC3 = 0; fNVtxSMRD = 0; fNVtxRestOfOffAxis = 0; fNVtxINGRID = 0; fNVtxOther = 0; fVertices->Clear(); const COMET::IHandle< COMET::IG4PrimaryVertexContainer >& vertices = event.Get("truth/G4PrimVertex00"); std::vector::const_iterator vertices_it = vertices->begin(); std::vector::const_iterator vertices_end; if(vertices->size() <= fMaxNVertices) { vertices_end = vertices->end(); } else { COMETWarn("Event contains " << vertices->size() << " vertices." << " Only the first " << fMaxNVertices << " will be saved."); vertices_end = vertices->begin(); vertices_end += fMaxNVertices; } for( ; vertices_it != vertices_end; ++vertices_it) { ITruthVertex* vtxToFill = new((*fVertices)[fNVtx]) ITruthVertex(); ++fNVtx; vtxToFill->ID = vertices_it->GetInteractionNumber(); vtxToFill->Position = vertices_it->GetPosition(); vtxToFill->Generator = vertices_it->GetGeneratorName(); vtxToFill->ReactionCode = vertices_it->GetReaction(); TLorentzVector pos = vertices_it->GetPosition(); gGeoManager->FindNode(pos.X(), pos.Y(), pos.Z()); std::string path = gGeoManager->GetPath(); vtxToFill->Subdetector = COMET::IoaAnalysisUtils::PathToSubdetector(path); // Fill vertex's fiducial distance and increment the relevant // vertex counter switch(vtxToFill->Subdetector) { case COMET::IoaAnalysisUtils::kFGD1: ++fNVtxFGD1; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kFGD2: ++fNVtxFGD2; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kP0D: ++fNVtxP0D; vtxToFill->Fiducial = IGeomInfo::P0D().WaterFiducialDistance(pos.X(), pos.Y(), pos.Z()) * unit::mm; break; case COMET::IoaAnalysisUtils::kDsECal: ++fNVtxDsECal; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kBrlECal: ++fNVtxBrlECal; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kP0DECal: ++fNVtxP0DECal; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kTPC1: ++fNVtxTPC1; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kTPC2: ++fNVtxTPC2; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kTPC3: ++fNVtxTPC3; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kMRD: ++fNVtxSMRD; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kOffAxis: ++fNVtxRestOfOffAxis; vtxToFill->Fiducial = -1.0; break; case COMET::IoaAnalysisUtils::kINGRID: ++fNVtxINGRID; vtxToFill->Fiducial = -1.0; break; default: ++fNVtxOther; vtxToFill->Fiducial = -1.0; break; } // END SWITCH ON SUBDETECTOR /// /// Fill primary particle trajectory information. /// These are particles leaving the vertex. std::vector::const_iterator primaries_it = vertices_it->GetPrimaryParticles().begin(); std::vector::const_iterator primaries_end = vertices_it->GetPrimaryParticles().end(); vtxToFill->NPrimaryTraj = 0; for( ; primaries_it != primaries_end; ++primaries_it) { vtxToFill->PrimaryTrajIDs.push_back( primaries_it->GetTrackId() ); ++(vtxToFill->NPrimaryTraj); } // For the info particles, it is a bit hard to specify the // behaviour without knowing exactly what the generator is // giving us. // For now, just look at particles with Id == -1 (-2 means // stuff going on inside the nucleus). Look for neutrinos and // non-neutrinos, assuming the latter are targets. // Loop through all info vertices and their particles until both // neutrino and target have been seen, then stop. bool found_neutrino = false; bool found_target = false; std::vector::const_iterator info_it = vertices_it->GetInfoVertex().begin(); std::vector::const_iterator info_end = vertices_it->GetInfoVertex().end(); for( ; info_it != info_end; ++info_it) { std::vector::const_iterator particles_it = info_it->GetPrimaryParticles().begin(); std::vector::const_iterator particles_end = info_it->GetPrimaryParticles().end(); Int_t temp_pdg; for( ; particles_it != particles_end; ++particles_it) { if(particles_it->GetTrackId() == -1) // is an incoming particle (including target) { temp_pdg = particles_it->GetPDGCode(); switch(temp_pdg) { case 12: // fall through case 14: // fall through case 16: // fall through case -12: // fall through case -14: // fall through case -16: vtxToFill->NeutrinoPDG = temp_pdg; vtxToFill->NeutrinoMomentum = particles_it->GetMomentum(); found_neutrino = true; break; default: vtxToFill->TargetPDG = temp_pdg; vtxToFill->TargetMomentum = particles_it->GetMomentum(); found_target = true; break; } if(found_neutrino && found_target) break; } } if(found_neutrino && found_target) break; } // END LOOP OVER INFO VERTICES } // END LOOP OVER VERTICES fVertices->Sort(); return true; } // // Truth Vertex // ===================================================================== // COMET::ITruthVerticesModule::ITruthVertex::ITruthVertex() : ID(-1), Position(-999999.9, -999999.9, -999999.9, -999999.9), Generator("Not Set"), ReactionCode("Not Set"), Subdetector(COMET::IoaAnalysisUtils::kNSubdetectors), Fiducial(-999999.9), NPrimaryTraj(-1), NeutrinoPDG(-1), NeutrinoMomentum(-999999.9, -999999.9, -999999.9, -999999.9), TargetPDG(-1), TargetMomentum(-999999.9, -999999.9, -999999.9, -999999.9) {} Int_t COMET::ITruthVerticesModule::ITruthVertex::Compare(const TObject* obj) const { // Sorry about the casting, but I don't want to slow down sorting // by checking that everything compared with is of the correct // type. In any case, you get what you ask for if you try and // use this with anything other than another vertex. if( this->ID < static_cast(obj)->ID ) { return -1; } else if( this->ID > static_cast(obj)->ID ) { return 1; } else // they must be equal { return 0; } }