|
//------------------------------------------------------------------------- // // BrFfsTrackingModule // // NEW tracking module for the FFS. The module will setup the module that // combines the local tracks found front and behind the D2 magnet. // The module does not instigate // the local tracking of T1 and T2. This is the responsibility of the // framework this will also allow a reconstruction job to substitute // the combine track part by another module without any penalty in // code organisation. // // //------------------------------------------------------------------------- // // $Id: BrFfsTrackingModule.cxx,v 1.6 2002/04/09 02:08:35 ouerdane Exp $ // // $Author: ouerdane $ // $Date: 2002/04/09 02:08:35 $ // $Copyright: 2000 Brahms Collaboration // //------------------------------------------------------------------------- #ifndef ROOT_TString #include "TString.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef ROOT_TMath #include "TMath.h" #endif #ifndef BRAT_BrFfsTrackingModule #include "BrFfsTrackingModule.h" #endif #ifndef BRAT_BrVector3D #include "BrVector3D.h" #endif #ifndef BRAT_BrModuleMatchTrack #include "BrModuleMatchTrack.h" #endif #ifndef BRAT_BrDetectorVolume #include "BrDetectorVolume.h" #endif #ifndef BRAT_BrMagnetVolume #include "BrMagnetVolume.h" #endif #ifndef BRAT_BrDetectorTrack #include "BrDetectorTrack.h" #endif #ifndef BRAT_BrSpectrometerTracks #include "BrSpectrometerTracks.h" #endif #ifndef BRAT_BrGeometryDbManager #include "BrGeometryDbManager.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrUnits #include "BrUnits.h" #endif #ifndef BRAT_BrRunInfoManager #include "BrRunInfoManager.h" #endif //__________________________________________________________________ ClassImp(BrFfsTrackingModule); //__________________________________________________________________ BrFfsTrackingModule::BrFfsTrackingModule() { //// default constructor DO NOT USE fCombineModule = 0; fTracksFrontBack = 0; SetNtuple(); // default is false } //__________________________________________________________________ BrFfsTrackingModule::BrFfsTrackingModule(const Char_t* Name, const Char_t* Title) : BrModule(Name, Title) { // ------------------------------------ // Normal constructor. // This Name & title is informational. // ------------------------------------ fCombineModule = new BrModuleMatchTrack("FFS", "FFS matching", "T1", "T2", "D2"); fTracksFrontBack = 0; SetNtuple(); // default is false fUseDb=kTRUE; } //__________________________________________________________________ BrFfsTrackingModule::~BrFfsTrackingModule() { //// default destructor if (fCombineModule) delete fCombineModule; } //__________________________________________________________________ void BrFfsTrackingModule::Init() { // ----------------------------------------- // Initialze the matching module and setup // local information for volumes // ----------------------------------------- fCombineModule->Init(); fCombineModule->SetDebugLevel(DebugLevel()); fCombineModule->SetVerbose(Verbose()); fFrontVolume = fCombineModule->GetFrontVolume(); fBackVolume = fCombineModule->GetBackVolume(); fMagnetVolume = fCombineModule->GetMagnetVolume(); // geometry BrGeometryDbManager *geoMan = BrGeometryDbManager::Instance(); fD1Volume = (BrMagnetVolume*)geoMan->GetDetectorVolume("BrMagnetVolume","D1"); fH1Volume = (BrDetectorVolume*)geoMan->GetDetectorVolume("BrDetectorVolume","TOF1"); // parameters BrParameterDbManager *parMan = BrParameterDbManager::Instance(); fH1Params = (BrDetectorParamsTof*)parMan-> GetDetectorParameters("BrDetectorParamsTof", "TOF1"); } //___________________________________________________________________ void BrFfsTrackingModule::SetMagnetFields() { // ---------------------------------------------------- // get field information from database // update info each time a field is needed (event loop) // ---------------------------------------------------- if(!fUseDb) return; BrRunInfoManager* runMan = BrRunInfoManager::Instance(); if (!runMan) { Abort("SetMagnetFields", "Could not get an instance of BrRunInfoManager"); return; } const BrRunInfo* run = runMan->GetCurrentRun(); if (run->GetRunNo() == -1) { Abort("SetMagnetFields", "Run number is -1, check it out"); return; } // get magnet settings fMagnetVolume->SetCurrent(run->GetD2Set(), run->GetD2Pol()); fD1Volume->SetCurrent(run->GetD1Set(), run->GetD1Pol()); if (Verbose() > 10) cout << " D1 field: " << fD1Volume->GetField() << " D2 field: " << fMagnetVolume->GetField() << endl; } //__________________________________________________________________________ void BrFfsTrackingModule::DefineHistograms() { // ----------------------------------------------------------- // define histograms for modules owned by this 'super module' // ----------------------------------------------------------- TDirectory* saveDir = gDirectory; TDirectory* histDir = gDirectory->mkdir("Tracking_FFS"); histDir->cd(); fCombineModule->Book(); fPathLength = new TH1F("length","Track length from BB vertex", 205, 500, 1000); fPartialLength = new TH1F("partialLength","Path length from H1 to D1 entrance", 1500, 500, 800); fTargetX = new TH1F("targetX","FFS target image (x)", 500, -25, 25); fTargetY = new TH1F("targetY","FFS target image (y)", 500, -25, 25); fFIncAng = new TH2F("fIncAng","Angle of incidence of front track on H1", 200, -20, 20, 180, 0, 15); fBIncAng = new TH2F("bIncAng","Angle of incidence of back track on H1", 200, -20, 20, 180, 0, 15); fAngCorr = new TH2F("incAngCorr", "Correlation between angle of incidence on H1" " of back track and front track", 180, 0, 15, 180, 0, 15); fH1CorrX = new TH2F("h1CorrX", "Correlation between back and front track " "projections on H1 X axis", 300, -30, 30, 300, -30, 30); fH1CorrY = new TH2F("h1CorrY", "Correlation between back and front track " "projections on H1 Y axis", 200, -10, 10, 200, -10, 10); fH1ProjX = new TH1F("h1ProjX", "Back track projection on H1 X axis", 300, -30, 30); fH1ProjY = new TH1F("h1ProjY", "Back track projection on H1 Y axis", 150, -15, 15); fH1Slat = new TH1F("h1Slat", "TOF1 slat pointed by back track", 60, -9.5, 50.5); fT1CorrX = new TH2F("t1CorrX", "Correlation between back and front track " "projections on T1 middle plane X axis", 300, -30, 30, 300, -30, 30); fT1CorrY = new TH2F("t1CorrY", "Correlation between back and front track " "projections on T1 middle plane Y axis", 200, -10, 10, 200, -10, 10); fT1DiffX = new TH1F("t1DiffX", "Correlation between back and front track " "projections on T1 middle plane X axis", 200, -10, 10); fT1DiffY = new TH1F("t1DiffY", "Correlation between back and front track " "projections on T1 middle plane Y axis", 200, -10, 10); fSwimStatus = new TH1F("swimStatD1", "D1 Swim Back Status: 0 = not good, 1 = Ok", 6, -0.5, 1.5); fMomentum= new TH1F("mom","Momentum of accepted FFS tracks",400,-20.,20.); fTrkX = new TH1F("trkVX", "Track intersection X with BB vertex plane", 400, -25, 25); fTrkY = new TH1F("trkVY", "Track intersection Y with BB vertex plane", 400, -25, 25); fTrkXY = new TH2F("trkVXY","Track projection on BB vertex plane", 200, -20, 20, 200, -20, 20); fCloseX = new TH1F("closeX", "Track X closest point to beam line", 400, -25, 25); fCloseY = new TH1F("closeY", "Track Y closest point to beam line", 400, -25, 25); fCloseXY = new TH2F("closeXY","Track Y vs X closest point to beam line", 200, -20, 20, 200, -20, 20); fCloseZ = new TH1F("closeZ", "Track Z closest point to beam line", 300, -150, 150); fCloseDiff = new TH1F("closeDiff", "Track Z closest - BB vtx", 600, -150, 150); fCloseCorr = new TH2F("closeCorr", "BB vtx vs Track closest Z to beamline", 300, -150, 150, 300, -150, 150); fStatus = new TH1F("trkStatus", "Ghostbusting track status", 15, 0.5, 5.5); fTheta = new TH1F("theta", "Theta in deg.", 600, 0, 45); fPhi = new TH1F("phi", "Phi in deg.", 600, -30, 30); if (fNtuple) fTrackInfo = new TNtuple("trackTree", "FFS Track Tree", "targX:targY:vtxX:vtxY:closeX:closeY:closeZ:" "momentum:theta:phi:d1ExitX:d1ExitY:d1EntrX:" "d1EntrY:Z0:status:frIncAngH1:bkIncAngH1:" "frProjH1:bkProjH1:bbMeth:d1Swim"); fNTracksAll = new TH1F("allTrks", "Number of matched tracks, all status", 11, -0.5, 10.5); fNTracks = new TH1F("goodTrks", "Number of matched tracks, status = 1", 11, -0.5, 10.5); gDirectory = saveDir; } //_________________________________________________________________________ void BrFfsTrackingModule::Begin() { // make sure we get the right field in for each run SetMagnetFields(); } //_________________________________________________________________________ void BrFfsTrackingModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // ------------------------------------------------------------- // Normal Event member function for FFS tracking // Result: BrDataTable "FfsTracks" added to outputnode. // The elements of this table are BrFfsTrack. // // The module has the following steps // a) combine Front Backs tracks // b) Track combined tracks back to nominal vertex // c) Add tracks to outNode // ------------------------------------------------------------- fCombineModule->Clear(); // ------- Get primary vertex: for now only BB vertex BrBbVertex* primVtx = (BrBbVertex*)inNode->GetObject("BB Vertex"); if (!primVtx) primVtx = (BrBbVertex*)outNode->GetObject("BB Vertex"); if (!primVtx) { if (Verbose() > 5) Warning("Event", "Could not find BB vertex"); // return; } // -------------------------------------------- // -- Combining local tracks -- // -------------------------------------------- BrEventNode* node = inNode; if (!inNode->GetDataTable("DetectorTrack T1") || !inNode->GetDataTable("DetectorTrack T2")) node = outNode; fCombineModule->CombineFrontBack(node); // ------------------------------------------------------ // -- Get matched tracks and complete tracking to vtx -- // ------------------------------------------------------ fTracksFrontBack = fCombineModule->GetMatchedTracks(); if(!fTracksFrontBack) { if (Verbose() > 5) Warning("Event", "No matched tracks array"); return; } if (!fTracksFrontBack->GetEntries()) return; // do tracking down to vertex BrDataTable *ffstable = new BrDataTable("FfsTracks"); outNode->AddDataTable(ffstable); TrackToVertex(ffstable, primVtx); // print out some stuff... for (Int_t i = 0; i < ffstable->GetEntries(); i++) { if (Verbose() > 10) ((BrFfsTrack*)ffstable->At(i))->Print(); else if (Verbose() > 15) ((BrFfsTrack*)ffstable->At(i))->Print("P"); else if (Verbose() > 20) ((BrFfsTrack*)ffstable->At(i))->Print("PD"); else if (Verbose() > 25) ((BrFfsTrack*)ffstable->At(i))->Print("A"); } } //___________________________________________________________________ void BrFfsTrackingModule::Finish() { // Finish entry method // Call combine module finish for statistics. fCombineModule->Finish(); } //___________________________________________________________________ void BrFfsTrackingModule::TrackToVertex(BrDataTable* ffstab, BrBbVertex* primVtx) { // ------------------------------------------- // Track back to vertex (here BB vertex) // ------------------------------------------- if(DebugLevel() > 0) cout << " In TrackToVertex(): " << " #Tracks = " << fTracksFrontBack->GetEntries() << endl; //Loop over tracks Int_t nTrk = fTracksFrontBack->GetEntries(); if (HistOn()) fNTracksAll->Fill(nTrk); Int_t nGood = 0; for (Int_t t = 0; t < nTrk; t++) { BrMatchedTrack *matchTrk = (BrMatchedTrack*)fTracksFrontBack->At(t); if (HistOn()) fStatus->Fill(matchTrk->GetStatus()); // check ghost tracks if (matchTrk->GetStatus() == 1) nGood++; // -------------- analyze best matched tracks // path length Double_t length = 0; //get the tracks in T1 & T2... BrDetectorTrack* t1Trk = matchTrk->GetFrontTrack(); BrDetectorTrack* t2Trk = matchTrk->GetBackTrack(); Double_t momentum = matchTrk->GetMomentum(); // Convert to global the T1 track line BrLine3D gTrkLine = fFrontVolume->LocalToGlobal(t1Trk->GetTrackLine()); // transform from global to D1 coordinates... BrLine3D d1ExTrkLine = fD1Volume->GlobalToLocal(gTrkLine); //project to exit of D1 (local to D1) BrPlane3D d1ExPlane = fD1Volume->GetEffectiveEdgeExitPlane(); BrVector3D d1ExPoint = d1ExPlane.GetIntersectionWithLine(d1ExTrkLine); d1ExTrkLine.SetOrigin(d1ExPoint); BrLine3D gD1ExTrkLine = fD1Volume->LocalToGlobal(d1ExTrkLine); BrVector3D gD1ExPoint = gD1ExTrkLine.GetOrigin(); // swim back track through D1 BrLine3D d1EnTrkLine = fD1Volume->SwimBack(d1ExTrkLine, momentum); // D1 swim back status // This is no good when momentum is 0.. // Need another method!! Bool_t swimStatus; if(TMath::Abs(momentum)<1000) swimStatus = fD1Volume->GetSwimStatus(d1EnTrkLine, d1ExTrkLine); else swimStatus=kFALSE; if (HistOn()) { if (swimStatus) fSwimStatus->Fill(1.); else fSwimStatus->Fill(0.); } // --------------------------------------------------------------- // convert entrance track line to global coord. sys. BrLine3D gD1EnTrkLine = fD1Volume->LocalToGlobal(d1EnTrkLine); Double_t vz = 0; if(primVtx) vz = primVtx->GetZ0(); // project trackIntoD1 to BB vertex BrVector3D vtx (0, 0, vz); BrVector3D norm(0, 0, 1); BrPlane3D vtxPlane(vtx, norm); BrVector3D trkVtx = vtxPlane.GetIntersectionWithLine(gD1EnTrkLine); // now, check the closest point between the beam line and the // track line BrLine3D beamLine(vtx, norm); BrLine3D closestLine(gD1EnTrkLine.GetShortestLineBetween(beamLine)); BrVector3D closePt = closestLine.GetOrigin(); if(HistOn()) { fTrkXY->Fill(trkVtx[0], trkVtx[1]); fTrkX->Fill(trkVtx[0]); fTrkY->Fill(trkVtx[1]); fCloseX->Fill(closePt[0]); fCloseY->Fill(closePt[1]); fCloseXY->Fill(closePt[0], closePt[1]); fCloseZ->Fill(closePt[2]); fCloseCorr->Fill(vz, closePt[2]); fCloseDiff->Fill(closePt[2] - vz); } // Calculate scattering angles and path length // FIX ME: assume here that the track comes from the BB vertex // I am not so sure about this (at least at this point // This should be reviews very soon. // BrVector3D vtx2d1 = gD1EnTrkLine.GetOrigin() - vtx; length += vtx2d1.Norm() + (gD1ExPoint - gD1EnTrkLine.GetOrigin()).Norm(); BrVector3D gD2EnPoint, gD2ExPoint; fMagnetVolume->LocalToGlobal(matchTrk->GetEntrance(), gD2EnPoint, 0); fMagnetVolume->LocalToGlobal(matchTrk->GetExit(), gD2ExPoint, 0); length += (gD2EnPoint - gD1ExPoint).Norm() + (gD2ExPoint - gD2EnPoint).Norm(); // scattering angles from bb vtx Float_t phi = (Float_t)TMath::ATan(vtx2d1[1] / vtx2d1[0]); Float_t xt = (Float_t)sqrt(vtx2d1[0]*vtx2d1[0] + vtx2d1[1]*vtx2d1[1]); Float_t theta = (Float_t)TMath::ATan(xt / vtx2d1[2] ); // ------- calculate track angle of incidence on tof plane // all here is in the TOF coord. system BrLine3D bkInc(BackTrackIncidence(t2Trk)); Double_t bkCosInc = bkInc.GetDirection().GetZ(); Double_t bkProjX = bkInc.GetOrigin().GetX(); Double_t bkProjY = bkInc.GetOrigin().GetY(); Double_t bkAngInc = TMath::ACos(bkCosInc)/BrUnits::degree; // line from d2 exit to H1 BrVector3D gH1Point; fH1Volume->LocalToGlobal(bkInc.GetOrigin(), gH1Point, 0); length += (gH1Point - gD2ExPoint).Norm(); // determine pointed slat and correct length for the slat pos Int_t slat = fH1Params->GetSlatNo(bkInc.GetOrigin()[0]); BrVector3D offset(0, 0, fH1Params->GetSlatPos(slat)[2]/bkCosInc); length += offset.Norm(); // diagnostic histograms if (HistOn()) { fH1Slat->Fill(slat); fPathLength->Fill(length); fH1ProjX->Fill(bkInc.GetOrigin()[0]); fH1ProjY->Fill(bkInc.GetOrigin()[1]); fTheta ->Fill(theta/BrUnits::degree); fPhi ->Fill(phi/BrUnits::degree); BrLine3D frInc(FrontTrackIncidence(t1Trk, momentum)); Double_t frCosInc = frInc.GetDirection().GetZ(); Double_t frProjX = frInc.GetOrigin().GetX(); Double_t frProjY = frInc.GetOrigin().GetY(); Double_t frAngInc = TMath::ACos(frCosInc)/BrUnits::degree; fBIncAng->Fill(bkProjX, bkAngInc); fFIncAng->Fill(frProjX, frAngInc); fAngCorr->Fill(frAngInc, bkAngInc); fH1CorrX->Fill(frProjX, bkProjX); fH1CorrY->Fill(frProjY, bkProjY); // check now where the back track line intersects T1 middle // plane BrVector3D bkInT1(BackTrackToT1(t2Trk, momentum)); BrVector3D frInT1 = t1Trk->GetPos(); fT1CorrX->Fill(frInT1.GetX(), bkInT1.GetX()); fT1CorrY->Fill(frInT1.GetY(), bkInT1.GetY()); fT1DiffX->Fill(bkInT1.GetX()-frInT1.GetX()); fT1DiffY->Fill(bkInT1.GetY()-frInT1.GetY()); // Assume, for now, Vertex is at 0,0,0 and define plane normal to // beam (i.e. Z = 0): BrPlane3D WrongVertex(1,0,0,0,1,0,0,0,0); BrVector3D IRIntersect = WrongVertex.GetIntersectionWithLine(gD1EnTrkLine); fTargetX->Fill(IRIntersect[0]); fTargetY->Fill(IRIntersect[1]); if (fNtuple) { Int_t s = 0; if (swimStatus) s = 1; // fill ntuple const Float_t ntVar[] = { IRIntersect[0], IRIntersect[1], trkVtx[0], trkVtx[1], closePt[0], closePt[1], closePt[2], momentum, theta, phi, d1ExPoint(0), d1ExPoint(1), gD1EnTrkLine.GetOrigin()[0], gD1EnTrkLine.GetOrigin()[1], vz , matchTrk->GetStatus(), frAngInc, bkAngInc, frProjX, bkProjX,1 , s }; fTrackInfo->Fill(ntVar); } } // create new FFS Track BrFfsTrack* trk = new BrFfsTrack(); trk->SetTrackId(t); trk->SetMatchedTrack(matchTrk); trk->SetD1EntranceLine(gD1EnTrkLine); trk->SetD1ExitLine(gD1ExTrkLine); trk->SetTheta(theta); trk->SetPhi(phi); trk->SetPointedSlat(slat); trk->SetProjOnTof(bkInc.GetOrigin() + offset); trk->SetPathLength(length); trk->SetTrackVertex(trkVtx); trk->SetD1SwimStatus(swimStatus); trk->SetD1H1PathLength(length - vtx2d1.Norm()); // fill table ffstab->Add(trk); if(HistOn()) { fMomentum->Fill(trk->GetMomentum()); fPartialLength->Fill(trk->GetPathLength() - vtx2d1.Norm()); } } if (HistOn()) fNTracks->Fill(nGood); } //___________________________________________________________________ BrLine3D BrFfsTrackingModule::BackTrackIncidence(BrDetectorTrack* track) { // ---------------------------------------------------------- // calculate projection on H1 plane and returns track line // with this proj. as the origin // ---------------------------------------------------------- BrPlane3D tofPlane(0,0,0, 0,1,0, 1,0,0); BrLine3D t2Line(track->GetTrackLine()); BrLine3D t2GLine(fBackVolume->LocalToGlobal(t2Line)); t2Line = fH1Volume->GlobalToLocal(t2GLine); BrVector3D proj = tofPlane.GetIntersectionWithLine(t2Line); t2Line.SetOrigin(proj); return t2Line; } //___________________________________________________________________ BrVector3D BrFfsTrackingModule::BackTrackToT1(BrDetectorTrack* track, Double_t P) { // ---------------------------------------------------------- // calculate projection on H1 plane and returns track line // with this proj. as the origin // ---------------------------------------------------------- BrPlane3D t1Plane(0,0,0, 0,1,0, 1,0,0); BrLine3D t2Line(track->GetTrackLine()); BrLine3D t2GLine(fBackVolume->LocalToGlobal(t2Line)); t2Line = fMagnetVolume->GlobalToLocal(t2GLine); BrLine3D newT2Line(fMagnetVolume->SwimBack(t2Line, P)); t2GLine = fMagnetVolume->LocalToGlobal(newT2Line); t2Line = fFrontVolume->GlobalToLocal(t2GLine); return t1Plane.GetIntersectionWithLine(t2Line);; } //___________________________________________________________________ BrLine3D BrFfsTrackingModule::FrontTrackIncidence(BrDetectorTrack* track, Double_t P) { // ---------------------------------------------------------- // calculate projection on H1 plane and returns track line // with this proj. as the origin // ---------------------------------------------------------- BrPlane3D tofPlane(0,0,0, 0,1,0, 1,0,0); BrLine3D t1Line(track->GetTrackLine()); BrLine3D t1GLine(fFrontVolume->LocalToGlobal(t1Line)); t1Line = fMagnetVolume->GlobalToLocal(t1GLine); BrLine3D newT1Line(fMagnetVolume->SwimForward(t1Line, P)); t1GLine = fMagnetVolume->LocalToGlobal(newT1Line); t1Line = fH1Volume->GlobalToLocal(t1GLine); BrVector3D proj = tofPlane.GetIntersectionWithLine(t1Line); t1Line.SetOrigin(proj); return t1Line; } //__________________________________________________________________________ // // $Log: BrFfsTrackingModule.cxx,v $ // Revision 1.6 2002/04/09 02:08:35 ouerdane // extend D4 matched tracks to D3 matched if no D3 matched tracks available // // Revision 1.5 2002/03/24 20:38:32 videbaek // Add finish entry // // Revision 1.4 2001/12/13 12:43:26 ouerdane // fixed a bug in the unit of theta and phi, added corresponding histos // // Revision 1.3 2001/12/07 21:36:59 videbaek // Add UseDb method to bypass RunDb. For particular Geant use. // Bypass (in FFS) swimstatus if momentum undetermined, but D1 is on. This happens for // some zero filed runs. SwinStatus would give an error message per track if not bypassed. // complete analysis in ffs even if no Vtx from BB. Set vz to 0 in Ntuples, etc. // Needed and useful for Geant + pp running when this comes up. // // Revision 1.2 2001/11/13 16:00:33 ouerdane // accept all tracks whatever ghost status // // Revision 1.1 2001/11/05 06:44:47 ouerdane // added new classes replacing older spectrometer track modules // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|