|
//--------------------------------------------------------------- // // BrTofPidModule // // Get Tof rdo hits and global tracks from input node // Get tof-track match table to pick up the right combination // Evaluate a PID (very rudimentary for the moment, we have to // evaluate the m2 resolution in function of the momentum, tof // and path length resolution) // //////////////////////////////////////////////////////////////// // // $Id: BrTofPidModule.cxx,v 1.29 2002/08/30 16:22:33 hagel Exp $ // $Author: hagel $ // $Date: 2002/08/30 16:22:33 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef ROOT_TMath #include "TMath.h" #endif #ifndef ROOT_TRandom #include "TRandom.h" #endif #ifndef ROOT_TString #include "TString.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef BRAT_BrTofPidModule #include "BrTofPidModule.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrGeometryDbManager #include "BrGeometryDbManager.h" #endif #ifndef BRAT_BrRunInfoManager #include "BrRunInfoManager.h" #endif #ifndef BRAT_BrBbVertex #include "BrBbVertex.h" #endif #ifndef BRAT_BrDataTable #include "BrDataTable.h" #endif #ifndef BRAT_BrTofRdo #include "BrTofRdo.h" #endif #ifndef BRAT_BrTofPid #include "BrTofPid.h" #endif #ifndef BRAT_BrTofTrackMatch #include "BrTofTrackMatch.h" #endif #ifndef BRAT_BrTd1TrackMatch #include "BrTd1TrackMatch.h" #endif #ifndef BRAT_BrTMrsFTrackMatch #include "BrTMrsFTrackMatch.h" #endif #ifndef BRAT_BrSpectrometerTracks #include "BrSpectrometerTracks.h" #endif #ifndef BRAT_BrDetectorTrack #include "BrDetectorTrack.h" #endif #ifndef BRAT_BrUnits #include "BrUnits.h" #endif #ifndef BRAT_BrPid #include "BrPid.h" #endif #ifndef BRAT_BrRichPid #include "BrRichPid.h" #endif #ifndef BRAT_BrC1Pid #include "BrC1Pid.h" #endif #ifndef BRAT_BrTd1Rdo #include "BrTd1Rdo.h" #endif #ifndef BRAT_BrEvent #include "BrEvent.h" #endif #ifndef BRAT_BrEventHeader #include "BrEventHeader.h" #endif #include <vector> //_____________________________________________________________ ClassImp(BrTofPidModule); //_____________________________________________________________ BrTofPidModule::BrTofPidModule() { // Default Constructor. (Do not use) fAccEvt = 0; SetDefaultParameters(); fInFfs = kFALSE; fInMrs = kFALSE; fInBfs = kFALSE; fMatchName = '0'; } //_____________________________________________________________ BrTofPidModule::BrTofPidModule(const Char_t *name, const Char_t *title) : BrModule(name, title) { fAccEvt = 0; SetDefaultParameters(); fInFfs = kFALSE; fInMrs = kFALSE; fInBfs = kFALSE; if (strcmp(name, "TOF1")==0) { fInFfs = kTRUE; fMatchName = "FFS Matching"; } else if (strcmp(name, "TOF2")==0) { fInBfs = kTRUE; fMatchName = "BFS Matching"; } else if (strcmp(name, "TOFW")==0) { fInMrs = kTRUE; fMatchName = "MRS Matching"; } else { Abort("BrTofPidModule", "Wrong detector name!! " "Should be TOF1 or TOF2 or TOFW"); return; } } //_____________________________________________________________ void BrTofPidModule::SetDefaultParameters() { // default parameters SetMaxDY(); // diff between track y pos and hit y pos (default 1cm) SetTrkOffsets(); // defaults are x = 0, y = 0, z = 0 SetTrkCuts(); // defaults are x = 10000, y = 10000, z = 10000 SetPScale(); // default is 1.00 SetNtuple(); // default is false SetUseBbVertex(); // default is true. SetUseTd1Time(); // default is false (pp only) SetUseTMrsTime(); // default is false (pp only) SetT0offset(); fFsLined = kFALSE; } //_____________________________________________________________ BrTofPidModule::~BrTofPidModule() { // destructor } //_____________________________________________________________ void BrTofPidModule::DefineHistograms() { // define histograms here TDirectory* saveDir = gDirectory; TDirectory* histDir = gDirectory->mkdir(Form("%sPid",GetName())); histDir->cd(); Int_t nSlats = fTofPar->GetNoSlats(); // histos if (fNtuple) { if (fInMrs) fPidTree = new TNtuple(Form("%sPid", GetName()), Form("%s Pid Tree", GetName()), "Z0:T0:bbM:tof:tofC:beta:momentum:ghost:dE:slat:trkL:mass2:" "vtxX:vtxY:vtxZ:theta:phi:eta:rap:pt:nPart:" "trkX:trkY:trkZ:hitX:hitY:hitZ:tofPan:seqno:trig"); else if (fInFfs) fPidTree = new TNtuple(Form("%sPid", GetName()), Form("%s Pid Tree", GetName()), "Z0:T0:bbM:tof:tofC:beta:momentum:ghost:d1Swim:dE:slat:trkL:mass2:" "vtxX:vtxY:vtxZ:theta:phi:eta:rap:pt:nPart:" "trkX:trkY:trkZ:hitX:hitY:hitZ:seqno:c1BlobE:trig:td1slat"); else if (fInBfs) fPidTree = new TNtuple(Form("%sPid", GetName()), Form("%s Pid Tree", GetName()), "Z0:T0:bbM:tof:tofC:beta:momentum:dE:slat:trkL:mass2:" "d1Swim:vtxX:vtxY:vtxZ:theta:phi:eta:rap:pt:nPart:" "trkX:trkY:trkZ:hitX:hitY:hitZ:seqno:richMass:trig"); } Float_t tmin, tmax, pmax; if (fInMrs) { if(fUseTMrsTime){ tmin = 8; tmax = 88; pmax = 4; } else { tmin = 13; tmax = 93; pmax = 4; } } else if (fInFfs) { if(fUseTd1Time) {tmin = 19; tmax = 29; pmax = 12; } else {tmin = 28; tmax = 37; pmax = 9; } } else if (fInBfs) { tmin = 58; tmax = 67; pmax = 15; } fPCut[0]=0; if (fInFfs) { fPCut[1] = 2; fPCut[2] = 3; fPCut[3] = 4; fPCut[4] = 5; fPCut[5] = 6; } if (fInBfs) { fPCut[1] = 4; fPCut[2] = 5; fPCut[3] = 6; fPCut[4] = 7; fPCut[5] = 8; } if (fInMrs) { fPCut[1] = 0.5; fPCut[2] = 1; fPCut[3] = 2; fPCut[4] = 3; fPCut[5] = 4; } // ----------- directories histDir->mkdir("No_Cuts"); histDir->mkdir("With_Cuts"); // ----------- Rough pid histDir->cd("No_Cuts"); fRoughTvVtx = new TH2F("hTvVtx", "Time of flight vs BB vertex", 600, -150, 150, 200, tmin-5, tmax+5); fRoughTvDE = new TH2F("hTvDE", "Time of flight vs #DeltaE", 300, 0.5, 2.5, 200, tmin, tmax); fRoughPid = new TH2F("hPid", "Time of flight vs Momentum, all tracks", 300, -pmax, pmax, 1000, tmin, tmax); fRoughDEvP = new TH2F("hDEvP", "#DeltaE vs Momentum", 300, -pmax, pmax, 300, 0.5, 4.5); fRoughM2[0] = new TH1F("hMass2", "Mass^{2} all momentum", 700, -1, 6); fRoughM2vP = new TH2F("hMass2vP", "Mass squared vs momentum", 500, -pmax, pmax, 700, -1, 6); fRoughBeta = new TH1F("hBeta", "Particle velocity #beta", 500, 0.7, 1.1); fRoughPvBeta = new TH2F("hPvsBeta", "Momentum vs 1/#beta", 250, -pmax, pmax, 500, 0.8, 1.8); fRoughTof = new TH1F("hTof", Form("Time of flight at %f m", fDist/100.), 800, tmin, tmax); fRoughDTrkTof = new TH2F("hDTrkTof", "Track-Tof #DeltaY vs #DeltaX", 400, -10, 10, 400, -10, 10); fRoughDy = new TH2F("hDy", "Track-Tof #DeltaY vs Slat", nSlats, 0.5, nSlats+0.5, 400, -10, 10); fRoughPartStat = new TH1F("hStat", "Number of clean particles per event", 10, -0.5, 9.5); fRoughLen = new TH1F("hPathLength", "Track length", 1000, 0, 2000); fRoughMom = new TH1F("hMomentum", "Track momentum", 400, -20, 20); fRoughEta = new TH1F("hEta", "Pseudo-rapidity #eta", 250, -10, 10); fRoughAccVtx = new TH1F("hVertex", "Primary vertex (event with identified particles)", 400, -200, 200); fRoughPtY = new TH2F("hPtY", "Pt vs Y", 500, -4, 4, 500, 0, 5); // ----------- Clean pid histDir->cd("With_Cuts"); fCleanTvVtx = new TH2F("hTvVtx", "Time of flight vs BB vertex", 600, -150, 150, 200, tmin-5, tmax+5); fCleanTvDE = new TH2F("hTvDE", "Time of flight vs #DeltaE", 300, 0.5, 2.5, 200, tmin, tmax); fCleanPid = new TH2F("hPid", "Time of flight vs Momentum, all tracks", 300, -pmax, pmax, 1000, tmin, tmax); fCleanDEvP = new TH2F("hDEvP", "#DeltaE vs Momentum", 300, -pmax, pmax, 300, 0.5, 4.5); fCleanM2[0] = new TH1F("hMass2", "Mass^{2} all momentum", 700, -1, 6); fCleanM2vP = new TH2F("hMass2vP", "Mass squared vs momentum", 500, -pmax, pmax, 700, -1, 6); fCleanBeta = new TH1F("hBeta", "Particle velocity #beta", 500, 0.7, 1.1); fCleanPvBeta = new TH2F("hPvsBeta", "Momentum vs 1/#beta", 250, -pmax, pmax, 500, 0.8, 1.8); fCleanTof = new TH1F("hTof", Form("Time of flight at %f m", fDist/100.), 800, tmin, tmax); fCleanDTrkTof = new TH2F("hDTrkTof", "Track-Tof #DeltaY vs #DeltaX", 400, -10, 10, 400, -10, 10); fCleanDy = new TH2F("hDy", "Track-Tof #DeltaY vs Slat", nSlats, 0.5, nSlats+0.5, 400, -10, 10); fCleanPartStat = new TH1F("hStat", "Number of clean particles per event", 10, -0.5, 9.5); fCleanLen = new TH1F("hPathLength", "Track length", 1000, 0, 2000); fCleanMom = new TH1F("hMomentum", "Track momentum", 400, -20, 20); fCleanEta = new TH1F("hEta", "Pseudo-rapidity #eta", 250, -10, 10); fCleanAccVtx = new TH1F("hVertex", "Primary vertex (event with identified particles)", 400, -200, 200); fCleanPtY = new TH2F("hPtY", "Pt vs Y", 500, -4, 4, 500, 0, 5); // ----------- histos platform dependent if (fInMrs) { histDir->cd("No_Cuts"); fRoughYvZ = new TH2F("hYvZ", "Track projection to vertex plane", 500, -25, 25, 500, -25, 25); fRoughM2[1] = new TH1F("hMass2Plus", "Mass^{2} momentum > 0", 700, -1, 6); fRoughM2[2] = new TH1F("hMass2minus", "Mass^{2} momentum < 0", 700, -1, 6); histDir->cd("With_Cuts"); fCleanYvZ = new TH2F("hYvZ", "Track projection to vertex plane", 500, -25, 25, 500, -25, 25); fCleanM2[1] = new TH1F("hMass2Plus", "Mass^{2} momentum > 0", 700, -1, 6); fCleanM2[2] = new TH1F("hMass2minus", "Mass^{2} momentum < 0", 700, -1, 6); } else { histDir->cd("No_Cuts"); fRoughYvX = new TH2F("hYvX", "Track projection to vertex plane", 500, -25, 25, 500, -25, 25); histDir->cd("With_Cuts"); fCleanYvX = new TH2F("hYvX", "Track projection to vertex plane", 500, -25, 25, 500, -25, 25); for (Int_t i = 1; i < 6; i++) { histDir->cd("No_Cuts"); fRoughM2[i] = new TH1F(Form("hMass2P%d", i), Form("Mass^{2}, %f < |P| < %f GeV/c", fPCut[i-1],fPCut[i]), 700, -1, 6); histDir->cd("With_Cuts"); fCleanM2[i] = new TH1F(Form("hMass2P%d", i), Form("Mass^{2}, %f < |P| < %f GeV/c", fPCut[i-1],fPCut[i]), 700, -1, 6); } } gDirectory = saveDir; } //_____________________________________________________________ void BrTofPidModule::Init() { // // initialize module // if(DebugLevel() > 1) cout << "Entering Init of BrTofPidModule" << endl; fSeqNo = 0; if (fInBfs) CheckFsAngles(); // parameter manager BrParameterDbManager* parDb = BrParameterDbManager::Instance(); fTofPar = (BrDetectorParamsTof*)parDb-> GetDetectorParameters("BrDetectorParamsTof", GetName()); // Setup reference distance here (and not in Constructor) // if(fInFfs){ if(fUseTd1Time) fDist=680; else fDist = 900; } if(fInBfs){ if(fUseTd1Time) fDist=1580; else fDist = 1800; } if(fInMrs) { if(fUseTMrsTime) fDist = 380.; else fDist = 433.5; } if (!fInMrs) return; // need this because the track projection on the Tof panel is // in global coordinate BrGeometryDbManager* geoDb = BrGeometryDbManager::Instance(); fPanelVol = new BrDetectorVolume* [fTofPar->GetNoPanels()]; for (Int_t p = 0; p < fTofPar->GetNoPanels(); p++) fPanelVol[p] = (BrDetectorVolume*)geoDb-> GetDetectorVolume("BrDetectorVolume", Form("TFP%d",p+1)); } //____________________________________________________________________ void BrTofPidModule::CheckFsAngles() { if (!fInBfs) return; // check if FFS and BFS are aligned BrRunInfoManager* runMan = BrRunInfoManager::Instance(); const BrRunInfo* run = runMan->GetCurrentRun(); if (run->GetRunNo() == -1) { Abort("Begin", "Current run number is -1. Aborting..."); return; } if (Verbose()) cout << " ---> FFS at " << run->GetFFSAngle() << " deg. " << " and BFS at " << run->GetBFSAngle() << " deg. " << endl; if (run->GetFFSAngle() == run->GetBFSAngle()) { fFsLined = kTRUE; if (Verbose()) cout << " ---> FS lined up " << endl; } } //___________________________________________________________ void BrTofPidModule::Begin() { // increment sequence number fSeqNo++; // sequence number fNoCleanPart = 0; // average particle within cuts per sequence fNoRoughPart = 0; // same whatever cuts fNoEvent = 0; // number of events in sequence // check FFS and BFS angles if (fInBfs) CheckFsAngles(); } //_____________________________________________________________ void BrTofPidModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // ------------------------------------------ // Event method to be called once per event. // 1- get primary vertex // 2- get hits and tracks // 3- get tof-track matching table // 4- pick up right tracks and hits // 5- get a preliminary pid // 6- store result in tables of BrTofPid // 7- fill histos before and after cuts // ------------------------------------------ fNoEvent++; // ------- get primary vertex (for now BB) BrBbVertex* vertex = (BrBbVertex*)inNode->GetObject("BB Vertex"); if (!vertex) vertex = (BrBbVertex*)outNode->GetObject("BB Vertex"); if(fUseBbVertex) if (!vertex) { if (Verbose() > 10) Warning("Event", " No beam beam stuff"); return; } Float_t Z0 = 0; Float_t T0 = 0; Int_t meth = 3; if(fUseBbVertex){ Z0 = vertex->GetZ0(); T0 = vertex->GetTime0(); meth = vertex->GetMethod(); } // ---------- get the matched tracks and hits table // matching objects BrDataTable* match = inNode->GetDataTable(fMatchName.Data()); if (!match) match = outNode->GetDataTable(fMatchName.Data()); if (!match) { if (Verbose() > 10) Warning("Event", "No tracks matching %s hits",GetName()); return; } // ------ tof table BrTofRdo *hits = (BrTofRdo*)inNode->GetObject(Form("TofRdo %s", GetName())); if (!hits) hits = (BrTofRdo*)outNode->GetObject(Form("TofRdo %s", GetName())); if (!hits) { if (Verbose() > 10) Warning("Event", "No reconstructed tof hits in this event"); return; } // ------ get track table BrDataTable* tracks = 0; if (fInFfs) tracks = inNode->GetDataTable("FfsTracks"); if (fInBfs) tracks = inNode->GetDataTable("FsTracks"); if (fInMrs) tracks = inNode->GetDataTable("MrsTracks"); if (!tracks) { if (fInFfs) tracks = outNode->GetDataTable("FfsTracks"); if (fInBfs) tracks = outNode->GetDataTable("FsTracks"); if (fInMrs) tracks = outNode->GetDataTable("MrsTracks"); } if (!tracks) { if (Verbose() > 10) Warning("Event", "No tracks at all in this event"); return; } BrDataTable* td1tab = 0; td1tab = inNode->GetDataTable("FfsMatchTd1"); BrDataTable* tmrstab = 0; tmrstab = inNode->GetDataTable("MrsFMatch"); // ------ get Cherenkov pid (only for histogramming) BrDataTable* chktab = 0; if (HistOn()) { if (fInFfs) chktab = inNode->GetDataTable("C1Pid"); if (fInBfs) chktab = inNode->GetDataTable("RichPid"); if (!chktab) { if (fInFfs) chktab = outNode->GetDataTable("C1Pid"); if (fInBfs) chktab = outNode->GetDataTable("RichPid"); } if (!fInMrs) if (!chktab) if (Verbose() > 25) Warning("Event", "No Cherenkov pid"); } // Find the start time in case TD! or TMrsF start counter are used // pp running // Int_t ntd1 = 0; if (fUseTd1Time) { if (td1tab) ntd1 = td1tab->GetEntries(); if (Verbose() > 10) for(Int_t k = 0; k < ntd1; k++) { BrTd1TrackMatch* td1m = (BrTd1TrackMatch*)td1tab->At(k); td1m->Print(); } } Int_t ntmrs = 0; if (fUseTMrsTime) { if (tmrstab) ntmrs = tmrstab->GetEntries(); if (Verbose() > 10) for(Int_t k = 0; k < ntmrs; k++) { BrTMrsFTrackMatch* tmrsm = (BrTMrsFTrackMatch*)tmrstab->At(k); tmrsm->Print(); } } if (Verbose() > 10) { cout << " -------------------------- " << endl << " Number of tracks : " << tracks->GetEntries() << endl << " Number of tof hits : " << hits->GetEntries() << endl << " Number of matching pairs : " << match->GetEntries() << endl; if (chktab) cout << " Number of cherenkov pid : " << chktab->GetEntries() << endl; if(tmrstab) cout << " Number of TMrs Matched : " << tmrstab->GetEntries() << endl; } // -------------------------------------------------------- // pick up right hit and track // we are now sure that the hit will correspond to a valid track Int_t nhit = hits->GetEntries(); // number of tof hits in event Int_t ntrk = tracks->GetEntries(); // number of tracks in event Int_t nchk = 0; if (chktab) nchk = chktab->GetEntries(); // number of chk pid in event Int_t nmch = match->GetEntries(); // number matching pairs if (HistOn()) fRoughPartStat->Fill(nmch); vector <BrTofRdo::BrTofHit*> tofh(nmch); // array of good hits vector <TObject*> track(nmch); // array of good tracks vector <BrPid*> chkpid(nmch); // chk pid object vector <BrTd1TrackMatch*> td1hit(nmch); // array of possible Td1Hits vector <BrTMrsFTrackMatch*> tmrshit(nmch); // array of possible TMrsFhits for (Int_t m = 0; m < nmch; m++) { Info(40,"Event","Starting loop, nhit = %d, ntrk = %d, nchk = %d, ntd1 = %dn",nhit,ntrk,nchk,ntd1); tofh [m] = 0; track [m] = 0; chkpid[m] = 0; td1hit[m] = 0; tmrshit[m] = 0; BrTofTrackMatch* pair = (BrTofTrackMatch*)match->At(m); // now select right hit thanks to the match object for (Int_t h = 0; h < nhit; h++) { BrTofRdo::BrTofHit* hit = hits->GetHit(h); if (hit->GetSlatno() != pair->GetHitId()) continue; tofh[m] = hit; } // now select right track proj thanks to the match object for (Int_t t = 0; t < ntrk; t++) { TObject* trk = tracks->At(t); if (trk->GetUniqueID() != pair->GetTrackId()) continue; track[m] = trk; } // now select chk pid objects if (chktab) for (Int_t p = 0; p < nchk; p++) { BrPid* chk = (BrPid*)chktab->At(p); if (chk->GetTrackId() != pair->GetTrackId()) continue; chkpid[m] = chk; } // // Lookup the Td1Hit corresponding to this track // if(fUseTd1Time){ if (td1tab){ Info(10,"Event","ntd1 = %d",ntd1); for (Int_t p = ntd1-1; p >=0; p--) { BrTd1TrackMatch* chk = (BrTd1TrackMatch*)td1tab->At(p); if (chk->GetTrackId() != pair->GetTrackId()) continue; td1hit[m] = chk; } } } // // Lookup the TMrsFHit corresponding to this track // if(fUseTMrsTime){ Info(10,"Event","ntmrs = %d",ntmrs); if (tmrstab){ for (Int_t p = ntmrs-1; p >=0; p--) { BrTMrsFTrackMatch* chk = (BrTMrsFTrackMatch*)tmrstab->At(p); if (chk->GetTrackId() != pair->GetTrackId()) continue; tmrshit[m] = chk; } } } } // vertex with track-tof trigger if (HistOn()) fRoughAccVtx->Fill(Z0); // check event trigger Int_t trig = CheckTrigger(inNode); // --------------------------------------------- // Particle identification Double_t c = BrUnits::c_light; // speed of light in cm/ns Int_t nCleanPart = 0; // number of particles // prepare output table BrDataTable* pidTab = new BrDataTable(Form("TofPid %s", GetName())); // estimate the number of particles that will be accepted for identification for (Int_t m = 0; m < nmch; m++) { if (tofh[m] == 0 || track[m] == 0) // in principle, not needed... continue; Float_t chksig = -2; if (chkpid[m]) { if (fInFfs) chksig = ((BrC1Pid*)chkpid[m])->GetBlobEnergy(); if (fInBfs) chksig = ((BrRichPid*)chkpid[m])->GetMass(); } // ------- get particle information Double_t length; Double_t mom; BrVector3D vtx; BrVector3D proj; Int_t ghost = 1; Int_t d1Swim = 1; Float_t theta; Float_t phi; if (fInMrs) { length = ((BrMrsTrack*)track[m])->GetPathLength(); mom = ((BrMrsTrack*)track[m])->GetMomentum(); ghost = ((BrMrsTrack*)track[m])->GetStatus(); vtx = ((BrMrsTrack*)track[m])->GetTrackVertex(); theta = ((BrMrsTrack*)track[m])->GetTheta(); phi = ((BrMrsTrack*)track[m])->GetPhi(); } if (fInFfs) { length = ((BrFfsTrack*)track[m])->GetPathLength(); mom = ((BrFfsTrack*)track[m])->GetMomentum(); vtx = ((BrFfsTrack*)track[m])->GetTrackVertex(); ghost = ((BrFfsTrack*)track[m])->GetStatus(); theta = ((BrFfsTrack*)track[m])->GetTheta(); phi = ((BrFfsTrack*)track[m])->GetPhi(); d1Swim = (Int_t)((BrFfsTrack*)track[m])->GetD1SwimStatus(); proj = ((BrFfsTrack*)track[m])->GetProjOnTof(); } if (fInBfs) { length = ((BrFsTrack*)track[m])->GetPathLength(); mom = ((BrFsTrack*)track[m])->GetMomentum(); vtx = ((BrFsTrack*)track[m])->GetTrackVertex(); theta = ((BrFsTrack*)track[m])->GetTheta(); phi = ((BrFsTrack*)track[m])->GetPhi(); d1Swim = (Int_t)((BrFsTrack*)track[m])->GetD1SwimStatus(); proj = ((BrFsTrack*)track[m])->GetProjOnTof2(); } Double_t tof = tofh[m]->GetTof(); Double_t ffsTrackLength = 0; Double_t bfsTrackLength = 0; Double_t td1TrackLength = 0; Float_t td1slat = 0; if (fUseTd1Time) { if (td1hit[m]) { T0 = td1hit[m]->GetTime(); tof = tofh[m]->GetRawTof() - T0 +fT0offset; if(fInFfs) { length = ((BrFfsTrack*)track[m])->GetPartialPath() + td1hit[m]->GetLength(); td1slat = td1hit[m]->GetSlat(); //Diagnostics (mainly used for bfs, but put here for symmetry) ffsTrackLength = ((BrFfsTrack*)track[m])->GetPartialPath(); td1TrackLength = td1hit[m]->GetLength(); } else if(fInBfs) { BrFfsTrack *ffsTrack = ((BrFsTrack*)track[m])->GetFfsTrack(); BrBfsTrack *bfsTrack = ((BrFsTrack*)track[m])->GetBfsTrack(); if(Verbose() > 30) { ((BrFsTrack*)track[m])->Dump(); if(ffsTrack) ffsTrack->Dump(); if(bfsTrack) bfsTrack->Dump(); } //We have to do following this way because there are some events //in which there is no ffsTrack in the BrFsTrack. I don't understand //that. Needs investigation. KH 25-APR-2002. In any case, //don't set td1slat if either track is missing. if(ffsTrack) ffsTrackLength = ffsTrack->GetPartialPath(); if(bfsTrack) bfsTrackLength = bfsTrack->GetPathLength(); td1TrackLength = td1hit[m]->GetLength(); length = ffsTrackLength + bfsTrackLength + td1TrackLength; if(ffsTrack && bfsTrack) td1slat = td1hit[m]->GetSlat(); Info(20,"Event","BFS track length = %f, td1slat = %d, ffsTrack = %d, bfsTrack = %d",length,td1slat,Int_t(ffsTrack),Int_t(bfsTrack)); } } } Float_t tmrsslat = 0; if (fUseTMrsTime) if (tmrshit[m]) { T0 = tmrshit[m]->GetTime(); tof = tofh[m]->GetRawTof() - T0 + fT0offset; length = ((BrMrsTrack*)track[m])->GetPartialPath() + tmrshit[m]->GetLength(); tmrsslat = 1; Info(20,"Event","Tmrshit = something, T0 = %f, tof = %f, length = %f",T0,tof,length); } Double_t ener = tofh[m]->GetEnergy(); Double_t beta = length/tof/c; Double_t tofC = fDist/beta/c; Double_t mass2 = mom*mom*(1/beta/beta - 1); if(fInMrs) { if ((fUseTMrsTime && tmrsslat > 0) || (!fUseTMrsTime)) { // ------ fill pid tab without any cleaning BrTofPid* pid = new BrTofPid(); pid->SetMatching((BrTofTrackMatch*)match->At(m)); pid->SetBeta(beta); pid->SetMomentum(mom); pid->SetDE(ener); pid->SetMass2(mass2); pidTab->Add(pid); Info(15,"Event","slat = %d, tof = %f, length = %f, p = %f, beta = %f, raw tof = %f, T0 = %f",tofh[m]->GetSlatNo(),tof,length,mom,beta,tofh[m]->GetRawTof(),T0); } } else if ((fUseTd1Time && td1slat > 0) || (!fUseTd1Time)) { // ------ fill pid tab without any cleaning BrTofPid* pid = new BrTofPid(); pid->SetMatching((BrTofTrackMatch*)match->At(m)); pid->SetBeta(beta); pid->SetMomentum(mom); pid->SetDE(ener); pid->SetMass2(mass2); pidTab->Add(pid); Info(15,"Event","slat = %d, tof = %f, length = %f, p = %f, beta = %f, raw tof = %f, T0 = %f",tofh[m]->GetSlatNo(),tof,length,mom,beta,tofh[m]->GetRawTof(),T0); } if (!HistOn()) continue; // --------------------------------------------- // histogram section // --------------------------------------------- // fill PID histo with everything (secondaries, bad tracks, etc) // pick up all info for histos BrTofTrackMatch* pair = (BrTofTrackMatch*)match->At(m); Float_t x = (vtx(0) - fTrkOffset[0])/fTrkCut[0]; Float_t y = (vtx(1) - fTrkOffset[1])/fTrkCut[1]; Float_t z = (vtx(2) - fTrkOffset[2] - Z0)/fTrkCut[2]; Float_t dist = x*x + y*y + z*z; Int_t slat = tofh[m]->GetSlatNum(); Int_t panel = pair->GetPanelId(); BrVector3D hitp = tofh [m]->GetLocalPos(); if (fInMrs) fPanelVol[panel]->GlobalToLocal(((BrMrsTrack*)track[m])->GetProjOnTof(), proj, 0); Float_t dX = proj(0) - hitp(0); Float_t dY = proj(1) - hitp(1); Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); Float_t rap = 0.5*TMath::Log((1+beta*TMath::Cos(theta))/(1-beta*TMath::Cos(theta))); Float_t pt = TMath::Abs(mom*TMath::Sin(theta)); // ----- fill ntuple with all the stuff if (fNtuple) { if (fInMrs) { const Float_t ntVar[] = { Z0, T0, meth, tof, tofC, beta, mom, ghost, ener, slat, length, mass2, vtx(0), vtx(1), vtx(2), theta, phi, eta, rap, pt, nmch, proj(0), proj(1), proj(2), hitp(0), hitp(1), hitp(2), panel, fSeqNo, trig }; fPidTree->Fill(ntVar); } else if (fInFfs) { const Float_t ntVar[] = { Z0, T0, meth, tof, tofC, beta, mom, ghost, d1Swim, ener, slat, length, mass2, vtx(0), vtx(1), vtx(2), theta, phi, eta, rap, pt, nmch, proj(0), proj(1), proj(2), hitp(0), hitp(1), hitp(2), fSeqNo, chksig, trig, td1slat }; fPidTree->Fill(ntVar); } else if (fInBfs) { const Float_t ntVar[] = { Z0, T0, meth, tof, tofC, beta, mom, ener, slat, length, mass2, d1Swim, vtx(0), vtx(1), vtx(2), theta, phi, eta, rap, pt, nmch, proj(0), proj(1), proj(2), hitp(0), hitp(1), hitp(2), fSeqNo, chksig, trig }; fPidTree->Fill(ntVar); } } // fill rough pid histos fRoughTvVtx ->Fill(Z0, tof); fRoughTvDE ->Fill(ener, tofC); fRoughPid ->Fill(mom, tofC); fRoughDEvP ->Fill(mom, ener); fRoughM2[0] ->Fill(mass2); fRoughM2vP ->Fill(mom, mass2); fRoughBeta ->Fill(beta); fRoughPvBeta ->Fill(mom, 1/beta); fRoughTof ->Fill(tofC); fRoughDTrkTof ->Fill(dX, dY); fRoughDy ->Fill(slat, dY); fRoughLen ->Fill(length); fRoughMom ->Fill(mom); fRoughEta ->Fill(eta); fRoughPtY ->Fill(rap, pt); if (fInMrs) { fRoughYvZ ->Fill(vtx(2) - Z0, vtx(1)); if(mom > 0) fRoughM2[1] ->Fill(mass2); else fRoughM2[2] ->Fill(mass2); } else { fRoughYvX ->Fill(vtx(0), vtx(1)); for (Int_t i = 1; i < 6; i++) if (TMath::Abs(mom) <= fPCut[i] && TMath::Abs(mom)> fPCut[i-1]) fRoughM2[i]->Fill(mass2); } if (fUseTd1Time) { if (td1slat == 0 || ghost !=1 || d1Swim !=1) continue; } else if (fUseTMrsTime) { if (tmrsslat == 0 || ghost !=1) continue; } else { // ---------- Clean Pid if (meth != 2) continue; // select only small tube vertex if (!IsWithinCuts(dist, dY, ghost, d1Swim)) continue; } nCleanPart++; // fill clean pid histos fCleanTvVtx ->Fill(Z0, tof); fCleanTvDE ->Fill(ener, tofC); fCleanPid ->Fill(mom, tofC); fCleanDEvP ->Fill(mom, ener); fCleanM2[0] ->Fill(mass2); fCleanM2vP ->Fill(mom, mass2); fCleanBeta ->Fill(beta); fCleanPvBeta ->Fill(mom, 1/beta); fCleanTof ->Fill(tofC); fCleanDTrkTof ->Fill(dX, dY); fCleanDy ->Fill(slat, dY); fCleanLen ->Fill(length); fCleanMom ->Fill(mom); fCleanEta ->Fill(eta); fCleanPtY ->Fill(rap, pt); if (fInMrs) { fCleanYvZ ->Fill(vtx(2) - Z0, vtx(1)); if(mom > 0) fCleanM2[1] ->Fill(mass2); else fCleanM2[2] ->Fill(mass2); } else { fCleanYvX ->Fill(vtx(0), vtx(1)); for (Int_t i = 1; i < 6; i++) if (TMath::Abs(mom) <= fPCut[i] && TMath::Abs(mom)>fPCut[i-1]) fCleanM2[i]->Fill(mass2); } } // -------------------------------- outNode->AddDataTable(pidTab); // accepted events (with at least one particle) fAccEvt++; fNoCleanPart += nCleanPart; fNoRoughPart += nmch; if (HistOn() && nCleanPart) { fCleanAccVtx->Fill(Z0); fCleanPartStat->Fill(nCleanPart); } } //_____________________________________________________________ Int_t BrTofPidModule::CheckTrigger(BrEventNode* node) { // private // check the trigger number // Have you checked what this does for overlap triggers ? // It just returns the HIGHEST trigger set. I do not think this is what you want. // Int_t trig = 0; if (node->IsA() != BrEvent::Class()) { Debug(25, "CheckTrigger", "Input node not a BrEvent"); return trig; } // Reinterpret as a BrEvent BrEvent* e = (BrEvent*)node; if (!e->IsEventHeader()) { Debug(25, "CheckTrigger", "No event header"); return trig; } BrEventHeader* header = e->GetEventHeader(); Int_t trigger = header->TriggerWord(1); trigger = trigger >> 8; for (Int_t i = 0; i < 8; i++) if (TESTBIT(trigger, i)) trig = i+1; Debug(25, "CheckTrigger", " Trigger is %d", trig); return trig; } //_____________________________________________________________ Bool_t BrTofPidModule::IsWithinCuts(Float_t vtxDist, Float_t dY, Int_t ghost, Int_t d1swim) { // private // check if particle candidate within users cuts if (vtxDist > 1 || TMath::Abs(dY) > fMaxDY || ghost != 1 || d1swim != 1) return kFALSE; return kTRUE; } //_____________________________________________________________ void BrTofPidModule::End() { // prints out some sequence statistics SetState(kEnd); if (Verbose() > 3) cout << " ---- " << GetName() << " PID - sequence # " << fSeqNo << "----" << endl << " Number of events analyzed : " << fNoEvent << endl << " Total number of particles without cuts : " << fNoRoughPart << endl << " Total number of particles within cuts : " << fNoCleanPart << endl << " Average particles per event (no cuts) : " << fNoRoughPart/fNoEvent << endl << " Average particles per event (within cuts) : " << fNoCleanPart/fNoEvent << endl << " --------------------------------------------------------" << endl; } //_____________________________________________________________ void BrTofPidModule::Finish() { SetState(kFinish); cout << " ---------------------------------------" << endl << " " << GetName() << " PID done." << endl << " Number of events with particle(s): " << fAccEvt << endl << " ---------------------------------------" << endl; } //_____________________________________________________________ void BrTofPidModule::Print(Option_t* option) const { // Standard information printout. // // Options: See BrModule::Print // TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) cout << endl << " Particle Identification Module" << endl << " Original author: Djamel Ouerdane" << endl << " Revisted by: $Author: hagel $" << endl << " Revision date: $Date: 2002/08/30 16:22:33 $" << endl << " Revision number: $Revision: 1.29 $ " << endl << endl << "*************************************************" << endl; } //_____________________________________________________________ // // // $Log: BrTofPidModule.cxx,v $ // Revision 1.29 2002/08/30 16:22:33 hagel // Add more verbose diagnostics // // Revision 1.28 2002/08/06 19:07:11 hagel // Format some parts according to ROOT coding convention (at least the part I like) // // Revision 1.27 2002/06/14 18:05:08 videbaek // Fix float int in strings and histograms names. // // Revision 1.26 2002/06/13 20:00:12 videbaek // Fix a typo. // Change fPcut to be floats so one can have p<1 cuts for Mrs. // // Revision 1.25 2002/06/13 18:54:31 videbaek // Add code to deal with TMrsF t0 startcounter for pp running. // Modify default length to 433.5 cm for MRS. Change some of summary plots // to have differential cuts on Mass2 rather than integral. // // Revision 1.24 2002/04/09 02:03:56 ouerdane // use BrFsTrack instead of BrBfsTrack for H2 PID // // Revision 1.23 2002/03/20 19:39:41 videbaek // Module now work with Td1 as start counter in addition to BB. // A seperate timing offset is needed for TD1, albeit if one creates // new slewings for pp I suspect this will be aligned properly to the // different patch length (680 cm). // // Revision 1.22 2002/01/31 17:10:40 ouerdane // added event trigger number in ntuple // // Revision 1.21 2001/12/20 17:41:49 ouerdane // introduced Cherenkov pid in Ntuples // // Revision 1.20 2001/12/13 14:25:17 ouerdane // changed the time range in the pid histo for MRS // // Revision 1.19 2001/12/13 13:38:39 ouerdane // cleaned module a lot, separeted histograms within or outside cuts, fill pid tables with every track-tof matches. The cuts are only for histograms // // Revision 1.18 2001/11/19 15:11:34 ouerdane // added all momentum information for H2 in the ntuple // // Revision 1.17 2001/11/12 18:51:35 ouerdane // corrected small bug in histos // // Revision 1.16 2001/11/07 10:32:35 ouerdane // updated matching module for H2 PID // // Revision 1.15 2001/11/05 07:05:20 ouerdane // changed a lot of things to deal with new track classes and BFS Pid // // Revision 1.14 2001/10/19 15:30:53 ouerdane // Added member fDist to fixe a path length for tof determination in PID plots // Fixed a couple of things in histograms and track length (added or substracted // a piece of length depending on the slat pos. in TOF1 // // Still a couple of things to be fixed but coming soon. // // Revision 1.13 2001/10/06 15:02:35 ouerdane // Fixed small error on TPM1 track vertex name // // Revision 1.12 2001/10/03 18:31:05 ouerdane // added more histograms for mass2 // // Revision 1.11 2001/10/03 17:20:08 ouerdane // minor error in DefineHistograms // // Revision 1.10 2001/10/03 17:19:19 ouerdane // added a scale factor for the momentum for testing // // Revision 1.9 2001/10/03 15:52:57 ouerdane // small error fixed on track cuts // // Revision 1.8 2001/10/02 20:19:33 ouerdane // remove member fNoPanels and setter method since it is now part of the // tof detector parameter class. // // Revision 1.7 2001/10/02 01:55:31 ouerdane // Updated the tof pid according to the tof track matching and primary vertex (cf email) // // Revision 1.6 2001/08/19 15:58:21 ouerdane // added event statistics check and vtx histogram for accepted events // // Revision 1.5 2001/07/24 14:41:10 ouerdane // Removed TList private member and added histograms as private members // (gain of time at run time) // Added a method public SetNtuple(Bool_t) , default is false // (an ntuple takes more and more space as the programs is running) // // Revision 1.4 2001/07/20 12:08:27 ouerdane // Changed + to - in diffX and diffY calculations when it comes // with the offsets between T2 and H1 in X and Y // // Revision 1.3 2001/06/29 13:51:53 cholm // Imported TOF classes from Djamel // // Revision 1.2 2001/06/21 13:01:14 ouerdane // Removed some obsolete things // added some table checks // // Revision 1.1 2001/06/19 12:53:21 ouerdane // Added class BrTofPidModule // It gets global tracks and reconstructed tof hits (BrTofRdo) // Match them and deduce a pid (still preliminary). // // This class was put in brat for the transition to brat 2 // but has not been tested in its present form (although it compiles // with brat). Don't use if you don't have a valid tof calibration // (cf. BrTofRdoModule) // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|