// ------------------------------------------------------------------------- // ----- CbmSttFitTracksQa source file ----- // ----- Created 30/03/06 by V. Friese ----- // ------------------------------------------------------------------------- #include "iostream.h" #include "TObjArray.h" #include "TVector3.h" #include "CbmMCTrack.h" #include "CbmRootManager.h" #include "CbmRuntimeDb.h" #include "CbmSttFitTracksQa.h" #include "CbmSttTrack.h" #include "CbmSttTrackMatch.h" #include "CbmSttGeomPoint.h" #include "CbmSttGeomHelix.h" #include "THelix.h" #include "TCanvas.h" #include "TView.h" #include "CbmSttHoughDefines.h" #include "CbmRunAna.h" #include "CbmField.h" #include "CbmSttVertex.h" #include "CbmSttHit.h" #define rangex 150 #define rangey 150 #define rangez 150 // ----- Default constructor ------------------------------------------- CbmSttFitTracksQa::CbmSttFitTracksQa() : CbmTask("STT Track Fitter QA", 1) { fEventNo = 0; CreateHistograms(); rootoutput = kFALSE; } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ CbmSttFitTracksQa::CbmSttFitTracksQa(Int_t iVerbose) : CbmTask("STTTrack Fitter QA", iVerbose) { fEventNo = 0; CreateHistograms(); if (iVerbose > 2) rootoutput = kTRUE; else rootoutput = kFALSE; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- CbmSttFitTracksQa::~CbmSttFitTracksQa() { } // ------------------------------------------------------------------------- // ----- Public method SetParContainers -------------------------------- void CbmSttFitTracksQa::SetParContainers() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- InitStatus CbmSttFitTracksQa::Init() { cout << "===========================================================" << endl;; cout << GetName() << ": Initialising..." << endl; // Get RootManager CbmRootManager* ioman = CbmRootManager::Instance(); if (! ioman) { cout << "-E- " << GetName() << "::Init: " << "RootManager not instantised!" << endl; return kFATAL; } // Get MCTrack array fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack"); if ( ! fMCTracks ) { cout << "-E- " << GetName() << "::Init: No MCTrack array!" << endl; return kFATAL; } // Get SttTrack array fSttTracks = (TClonesArray*) ioman->GetObject("STTTrack"); if ( ! fSttTracks ) { cout << "-E- " << GetName() << "::Init: No SttTrack array!" << endl; return kERROR; } // Get SttHit array fSttHits = (TClonesArray*) ioman->GetObject("STTHit"); if ( ! fSttHits ) { cout << "-E- " << GetName() << "::Init: No SttHit array!" << endl; return kERROR; } // Get SttTrackMatch array fMatches = (TClonesArray*) ioman->GetObject("STTTrackMatch"); if ( ! fMatches ) { cout << "-E- " << GetName() << "::Init: No SttTrackMatch array!" << endl; return kERROR; } // Get SttVertex array fSttVertices = (TClonesArray*) ioman->GetObject("STTVertex"); if ( ! fSttTracks ) { cout << "-E- " << GetName() << "::Init: No SttVertex array found!" << endl; return kERROR; } fTargetPos.SetXYZ(0., 0., 0.); // Output cout << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") " << endl; if (fActive) cout << " ***** Task is ACTIVE *****" << endl; cout << "===========================================================" << endl << endl; return kSUCCESS; } // ------------------------------------------------------------------------- // ----- Public method ReInit ------------------------------------------ InitStatus CbmSttFitTracksQa::ReInit() { cout << "===========================================================" << endl;; cout << GetName() << ": Reinitialising..." << endl; fTargetPos.SetXYZ(0., 0., 0.); // Output cout << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") " << endl; if (fActive) cout << " ***** Task is ACTIVE *****" << endl; cout << "===========================================================" << endl << endl; return kSUCCESS; } // ------------------------------------------------------------------------- // ----- Public method Exec -------------------------------------------- void CbmSttFitTracksQa::Exec(Option_t* opt) { fEventNo++; cout << "****************************** qa *************************" << endl; // Clear matching map and reset eventwise counters Int_t nRec = fSttTracks->GetEntriesFast(), nMtc = fMatches->GetEntriesFast(); if ( nMtc != nRec ) { cout << "-E- " << GetName() << ":Exec Number of SttMatches (" << nMtc << ") does not equal number of SttTracks (" << nRec << ")" << endl; cout << " Skipping event " << endl; return; } CbmSttGeomHelix mcProtonHelix, recoProtonHelix; CbmMCTrack *mcProton; CbmSttTrack *recoProton; TCanvas *qaCanvas; TView *view; if (rootoutput) { qaCanvas = new TCanvas("hel", "hel"); qaCanvas->Show(); qaCanvas->Update(); view = new TView(1); view->SetRange(-rangex, -rangey, -rangez, rangex, rangey, rangez); drawDetector(); drawAxis(); qaCanvas->Update(); } Int_t nRecAboveFourHits = 0; Bool_t eventHasAnnihilation = kFALSE, eventHasPionDecay = kFALSE, antiProtonFound = kFALSE, normalProtonFound = kFALSE, positivePionFound = kFALSE, negativePionFound = kFALSE; // loop over all reconstructed tracks for (Int_t iRec=0; iRec < nRec; iRec++) { cout << "irec: " << iRec << " nrec: " << nRec << endl; CbmSttTrackMatch *match = (CbmSttTrackMatch*) fMatches->At(iRec); CbmSttTrack *track = (CbmSttTrack*) fSttTracks->At(iRec); if ( ! match ) { cout << "-W- " << GetName() << "::Exec: Empty SttTrackMatch at " << iRec << endl; continue; } if ( ! track ) { cout << "-W- " << GetName() << "::Exec: Empty SttTrack at " << iRec << endl; continue; } Int_t iMC = match->GetMCTrackID(); if (iMC == -1 ) { cout << "-I- " << GetName() << ": SttTrack " << iRec << " has no corresponding MCTrack " << endl; continue; } CbmMCTrack *mcTrack = (CbmMCTrack*) fMCTracks->At(iMC); if ( ! mcTrack ) { cout << "-W- " << GetName() << "::Exec: Empty MCTrack at " << iMC << endl; continue; } CbmSttGeomHelix recoHelix(track), mcHelix(mcTrack); CbmSttGeomPoint vertex(mcTrack->GetStartVertex().X(), mcTrack->GetStartVertex().Y(), mcTrack->GetStartVertex().Z()); TVector3 momentum = recoHelix.GetMomentumAt(vertex); if (rootoutput) { recoHelix.DrawFromVertex(vertex, 0.5, 2); mcHelix.DrawFromVertex(vertex, 0.5, 3); drawHits(iRec); } fLambdaHitsPerTrack->Fill(track->GetNHits()); fLambdaParticleTypes->Fill(mcTrack->GetPdgCode()); if (track->GetNHits() > 4) { nRecAboveFourHits++; } switch (mcTrack->GetPdgCode()) { case -2212: fLambdaHitsPerTrackAntiproton->Fill(track->GetNHits()); break; case -211: fLambdaHitsPerTrackPiminus->Fill(track->GetNHits()); break; case -11: fLambdaHitsPerTrackEplus->Fill(track->GetNHits()); break; case -13: fLambdaHitsPerTrackMuplus->Fill(track->GetNHits()); break; case 2212: fLambdaHitsPerTrackProton->Fill(track->GetNHits()); break; case 211: fLambdaHitsPerTrackPiplus->Fill(track->GetNHits()); break; case 11: fLambdaHitsPerTrackEminus->Fill(track->GetNHits()); break; case 13: fLambdaHitsPerTrackMuminus->Fill(track->GetNHits()); break; default: fLambdaHitsPerTrackUnknown->Fill(track->GetNHits()); break; } CbmMCTrack *mcTrackMother; if (mcTrack->GetMotherID() != -1) { mcTrackMother = (CbmMCTrack*) fMCTracks->At(mcTrack->GetMotherID()); fLambdaParticleOrigins->Fill(mcTrackMother->GetPdgCode()); if (mcTrackMother->GetPdgCode() == -2212) { eventHasAnnihilation = kTRUE; } if (mcTrackMother->GetPdgCode() == -211 || mcTrackMother->GetPdgCode() == 211) { eventHasPionDecay = kTRUE; } switch (mcTrack->GetPdgCode()) { case -2212: fLambdaParticleOriginsAntiproton->Fill(mcTrackMother->GetPdgCode()); break; case -211: fLambdaParticleOriginsPiminus->Fill(mcTrackMother->GetPdgCode()); break; case -11: fLambdaParticleOriginsEplus->Fill(mcTrackMother->GetPdgCode()); break; case -13: fLambdaParticleOriginsMuplus->Fill(mcTrackMother->GetPdgCode()); break; case 2212: fLambdaParticleOriginsProton->Fill(mcTrackMother->GetPdgCode()); break; case 211: fLambdaParticleOriginsPiplus->Fill(mcTrackMother->GetPdgCode()); break; case 11: fLambdaParticleOriginsEminus->Fill(mcTrackMother->GetPdgCode()); break; case 13: fLambdaParticleOriginsMuminus->Fill(mcTrackMother->GetPdgCode()); break; default: fLambdaParticleOriginsUnknown->Fill(mcTrackMother->GetPdgCode()); break; } } else { fLambdaParticleOrigins->Fill(0); } if ((mcTrack->GetPdgCode() == 2212) && (track->GetNHits() > 4)) { mcProton = mcTrack; recoProton = track; mcProtonHelix = mcHelix; normalProtonFound = kTRUE; } } // Loop over all reconstructed tracks // the number of tracks reconstructed fLambdaRecoTracks->Fill(nRec); // number of reconstructed tracks in pbar annihilation events if (eventHasAnnihilation) fLambdaRecoTracksAnnihilation->Fill(nRecAboveFourHits); // number of reconstructed tracks in pion decay events if (eventHasPionDecay) fLambdaRecoTracksPionDecay->Fill(nRecAboveFourHits); // number of reconstructed tracks after selecting tracks with more than 4 hits fLambdaRecoTracksAbove->Fill(nRecAboveFourHits); // if all tracks reconstructed // if (antiProtonFound && normalProtonFound && positivePionFound && negativePionFound && nRecAboveFourHits == 4) if (normalProtonFound) { CbmSttGeomHelix recoProtonHelix(recoProton); // four hits selected fLambdaRecoTracksGood->Fill(nRecAboveFourHits); // if charge sign is reconstructed correctly if (recoProtonHelix.GetChargeSign() > 0.) { // chisquares per track fLambdaProtonChiSquareLong->Fill(recoProton->GetChi2Long()); fLambdaProtonChiSquareRad->Fill(recoProton->GetChi2Rad()); // number of vertices reconstructed //fLambdaRecoVertices->Fill(fSttVertices->GetEntriesFast()); // montecarlo distribution of vertices lambda Double_t mcVertexProtonX = mcProton->GetStartVertex().X(), mcVertexProtonY = mcProton->GetStartVertex().Y(), mcVertexProtonZ = mcProton->GetStartVertex().Z(); Double_t mcVertexProtonDistance = sqrt(mcVertexProtonX * mcVertexProtonX + mcVertexProtonY * mcVertexProtonY + mcVertexProtonZ * mcVertexProtonZ); fLambdaMCVertexProton->Fill(mcVertexProtonDistance); // decay particle momenta CbmSttGeomPoint origin(0., 0., 0.); TVector3 recoProtonMomentum = recoProtonHelix.GetMomentumAt(origin), mcProtonMomentum = mcProtonHelix.GetMomentumAt(origin); fLambdaProtonMomentumX->Fill(recoProtonMomentum.X() - mcProtonMomentum.X()); fLambdaProtonMomentumY->Fill(recoProtonMomentum.Y() - mcProtonMomentum.Y()); fLambdaProtonMomentumZ->Fill(recoProtonMomentum.Z() - mcProtonMomentum.Z()); fLambdaProtonMomentumPt->Fill(recoProtonHelix.GetPt() - mcProtonHelix.GetPt()); fLambdaProtonMomentumP->Fill(recoProtonHelix.GetP() - mcProtonHelix.GetP()); cout << "[[[[[[[[[[[[[reco: " << recoProtonHelix.GetPt() / 0.006 << " " << mcProtonHelix.GetPt() / 0.006 << endl; } // charge signs reconstructed okay } // all four tracks reconstructed if (rootoutput) { qaCanvas->Update(); char waitchar; cout << "Press any key to continue." << endl; cin >> waitchar; delete view; delete qaCanvas; } } // ------------------------------------------------------------------------- void CbmSttFitTracksQa::CreateHistograms() { fLambdaVertexMomentumAngleLambda = new TH1F("fLambdaVertexMomentumAngleLambda", "fLambdaVertexMomentumAngleLambda", 100, -0.5, 0.5); fLambdaVertexMomentumAngleAntiLambda = new TH1F("fLambdaVertexMomentumAngleAntiLambda", "fLambdaVertexMomentumAngleAntiLambda", 100, -0.5, 0.5); fLambdaBadPBarTheta = new TH1F("fLambdaBadPBarTheta", "fLambdaBadPBarTheta", 100, -2 * dPi, 2 * dPi); fLambdaGoodPBarTheta = new TH1F("fLambdaGoodPBarTheta", "fLambdaGoodPBarTheta", 100, -2 * dPi, 2 * dPi); fLambdaBadPBarHits = new TH1F("fLambdaBadPBarHits", "fLambdaBadPBarHits", 100, -0.5, 99.5); fLambdaGoodPBarHits = new TH1F("fLambdaGoodPBarHits", "fLambdaGoodPBarHits", 100, -0.5, 99.5); fLambdaSelectorHits = new TH1F("fLambdaSelectorHits", "fLambdaSelectorHits", 10, -0.5, 9.5); fLambdaSelectorHitsBad = new TH1F("fLambdaSelectorHitsBad", "fLambdaSelectorHitsBad", 10, -0.5, 9.5); fLambdaProtonChiSquareLong = new TH1F("fLambdaProtonChiSquareLong", "fLambdaProtonChiSquareLong", 40, 0.0, 100.); fLambdaPiPlusChiSquareLong = new TH1F("fLambdaPiPlusChiSquareLong", "fLambdaPiPlusChiSquareLong", 40, 0.0, 100.); fLambdaPiMinusChiSquareLong = new TH1F("fLambdaPiMinusChiSquareLong", "fLambdaPiMinusChiSquareLong", 40, 0.0, 100.); fLambdaAntiProtonChiSquareLong = new TH1F("fLambdaAntiProtonChiSquareLong", "fLambdaAntiProtonChiSquareLong", 40, 0.0, 100.); fLambdaProtonChiSquareRad = new TH1F("fLambdaProtonChiSquareRad", "fLambdaProtonChiSquareRad", 40, 0.0, 100.); fLambdaPiPlusChiSquareRad = new TH1F("fLambdaPiPlusChiSquareRad", "fLambdaPiPlusChiSquareRad", 40, 0.0, 100.); fLambdaPiMinusChiSquareRad = new TH1F("fLambdaPiMinusChiSquareRad", "fLambdaPiMinusChiSquareRad", 40, 0.0, 100.); fLambdaAntiProtonChiSquareRad = new TH1F("fLambdaAntiProtonChiSquareRad", "fLambdaAntiProtonChiSquareRad", 40, 0.0, 100.); fLambdaCosThetaLambdaMc = new TH1F("fLambdaCosThetaLambdaMc", "fLambdaCosThetaLambdaMc", 20, -1.0, 1.0); fLambdaCosThetaLambdaReco = new TH1F("fLambdaCosThetaLambdaReco", "fLambdaCosThetaLambdaReco", 20, -1.0, 1.0); fLambdaCosThetaAntiLambdaReco = new TH1F("fLambdaCosThetaAntiLambdaReco", "fLambdaCosThetaAntiLambdaReco", 20, -1.0, 1.0); fLambdaPhiLambdaReco = new TH1F("fLambdaPhiLambdaReco", "fLambdaPhiLambdaReco", 20, -dPi, dPi); fLambdaPhiAntiLambdaReco = new TH1F("fLambdaPhiAntiLambdaReco", "fLambdaPhiAntiLambdaReco", 20, -dPi, dPi); fLambdaOpeningPhiLambdaLambdabarReco = new TH1F("fLambdaOpeningPhiLambdaLambdabarReco", "fLambdaOpeningPhiLambdaLambdabarReco", 20, -2 * dPi, 2 * dPi); fLambdaOpeningThetaLambdaLambdabarReco = new TH1F("fLambdaOpeningThetaLambdaLambdabarReco", "fLambdaOpeningThetaLambdaLambdabarReco", 20, -2 * dPi, 2 * dPi); fLambdaTotalPxReco = new TH1F("fLambdaTotalPxReco", "fLambdaTotalPxReco", 300, -1.0, 1.0); fLambdaTotalPyReco = new TH1F("fLambdaTotalPyReco", "fLambdaTotalPyReco", 300, -1.0, 1.0); fLambdaTotalPzReco = new TH1F("fLambdaTotalPzReco", "fLambdaTotalPzReco", 300, 0.0, 4.0); fLambdaTotalEnReco = new TH1F("fLambdaTotalEnReco", "fLambdaTotalEnReco", 300, 0.0, 4.0); fLambdaTotalPxMc = new TH1F("fLambdaTotalPxMc", "fLambdaTotalPxMc", 300, -1.0, 1.0); fLambdaTotalPyMc = new TH1F("fLambdaTotalPyMc", "fLambdaTotalPyMc", 300, -1.0, 1.0); fLambdaTotalPzMc = new TH1F("fLambdaTotalPzMc", "fLambdaTotalPzMc", 300, 0.0, 4.0); fLambdaTotalEnMc = new TH1F("fLambdaTotalEnMc", "fLambdaTotalEnMc", 300, 0.0, 4.0); fLambdaInvMassLambda = new TH1F("fLambdaInvMassLambda", "fLambdaInvMassLambda", 1000, 0.0, 2.0); fLambdaInvMassAntiLambda = new TH1F("fLambdaInvMassAntiLambda", "fLambdaInvMassAntiLambda", 1000, 0.0, 2.0); fLambdaInvMassSum = new TH1F("fLambdaInvMassSum", "fLambdaInvMassSum", 1000, 0.0, 2.0); fLambdaLambdaMomentumX = new TH1F("fLambdaLambdaMomentumX", "fLambdaLambdaMomentumX", 100, -0.5, 0.5); fLambdaLambdaMomentumY = new TH1F("fLambdaLambdaMomentumY", "fLambdaLambdaMomentumY", 100, -0.5, 0.5); fLambdaLambdaMomentumZ = new TH1F("fLambdaLambdaMomentumZ", "fLambdaLambdaMomentumZ", 100, -0.5, 0.5); fLambdaAntiLambdaMomentumX = new TH1F("fLambdaAntiLambdaMomentumX", "fLambdaAntiLambdaMomentumX", 100, -0.5, 0.5); fLambdaAntiLambdaMomentumY = new TH1F("fLambdaAntiLambdaMomentumY", "fLambdaAntiLambdaMomentumY", 100, -0.5, 0.5); fLambdaAntiLambdaMomentumZ = new TH1F("fLambdaAntiLambdaMomentumZ", "fLambdaAntiLambdaMomentumZ", 100, -0.5, 0.5); fLambdaPiPlusMomentumX = new TH1F("fLambdaPiPlusMomentumX", "fLambdaPiPlusMomentumX", 100, -0.5, 0.5); fLambdaPiPlusMomentumY = new TH1F("fLambdaPiPlusMomentumY", "fLambdaPiPlusMomentumY", 100, -0.5, 0.5); fLambdaPiPlusMomentumZ = new TH1F("fLambdaPiPlusMomentumZ", "fLambdaPiPlusMomentumZ", 100, -0.5, 0.5); fLambdaPiPlusMomentumPt = new TH1F("fLambdaPiPlusMomentumPt", "fLambdaPiPlusMomentumPt", 100, -0.5, 0.5); fLambdaPiPlusMomentumP = new TH1F("fLambdaPiPlusMomentumP", "fLambdaPiPlusMomentumP", 100, -0.5, 0.5); fLambdaPiMinusMomentumX = new TH1F("fLambdaPiMinusMomentumX", "fLambdaPiMinusMomentumX", 100, -0.5, 0.5); fLambdaPiMinusMomentumY = new TH1F("fLambdaPiMinusMomentumY", "fLambdaPiMinusMomentumY", 100, -0.5, 0.5); fLambdaPiMinusMomentumZ = new TH1F("fLambdaPiMinusMomentumZ", "fLambdaPiMinusMomentumZ", 100, -0.5, 0.5); fLambdaPiMinusMomentumPt = new TH1F("fLambdaPiMinusMomentumPt", "fLambdaPiMinusMomentumPt", 100, -0.5, 0.5); fLambdaPiMinusMomentumP = new TH1F("fLambdaPiMinusMomentumP", "fLambdaPiMinusMomentumP", 100, -0.5, 0.5); fLambdaProtonMomentumX = new TH1F("fLambdaProtonMomentumX", "fLambdaProtonMomentumX", 100, -0.5, 0.5); fLambdaProtonMomentumY = new TH1F("fLambdaProtonMomentumY", "fLambdaProtonMomentumY", 100, -0.5, 0.5); fLambdaProtonMomentumZ = new TH1F("fLambdaProtonMomentumZ", "fLambdaProtonMomentumZ", 100, -0.5, 0.5); fLambdaProtonMomentumPt = new TH1F("fLambdaProtonMomentumPt", "fLambdaProtonMomentumPt", 100, -0.5, 0.5); fLambdaProtonMomentumP = new TH1F("fLambdaProtonMomentumP", "fLambdaProtonMomentumP", 100, -0.5, 0.5); fLambdaAntiProtonMomentumX = new TH1F("fLambdaAntiProtonMomentumX", "fLambdaAntiProtonMomentumX", 100, -0.5, 0.5); fLambdaAntiProtonMomentumY = new TH1F("fLambdaAntiProtonMomentumY", "fLambdaAntiProtonMomentumY", 100, -0.5, 0.5); fLambdaAntiProtonMomentumZ = new TH1F("fLambdaAntiProtonMomentumZ", "fLambdaAntiProtonMomentumZ", 100, -0.5, 0.5); fLambdaAntiProtonMomentumPt = new TH1F("fLambdaAntiProtonMomentumPt", "fLambdaAntiProtonMomentumPt", 100, -0.5, 0.5); fLambdaAntiProtonMomentumP = new TH1F("fLambdaAntiProtonMomentumP", "fLambdaAntiProtonMomentumP", 100, -0.5, 0.5); fLambdaVertexResolutionX = new TH1F("fLambdaVertexResolutionX", "fLambdaVertexResolutionX", 50, -10.0, 10.0); fLambdaVertexResolutionY = new TH1F("fLambdaVertexResolutionY", "fLambdaVertexResolutionY", 50, -10.0, 10.0); fLambdaVertexResolutionZ = new TH1F("fLambdaVertexResolutionZ", "fLambdaVertexResolutionZ", 50, -10.0, 10.0); fLambdaVertexResolution = new TH1F("fLambdaVertexResolution", "fLambdaVertexResolution", 50, -10.0, 10.0); fLambdaRecoTracks = new TH1F("fLambdaRecoTracks", "fLambdaRecoTracks", 20, -0.5, 19.5); fLambdaHitsPerTrack = new TH1F("fLambdaHitsPerTrack", "fLambdaHitsPerTrack", 100, -0.5, 99.5); fLambdaRecoTracksAbove = new TH1F("fLambdaRecoTracksAbove", "fLambdaRecoTracksAbove", 20, -0.5, 19.5); fLambdaRecoTracksAnnihilation = new TH1F("fLambdaRecoTracksAnnihilation", "fLambdaRecoTracksAnnihilation", 20, -0.5, 19.5); fLambdaRecoTracksPionDecay = new TH1F("fLambdaRecoTracksPionDecay", "fLambdaRecoTracksPionDecay", 20, -0.5, 19.5); fLambdaRecoTracksGood = new TH1F("fLambdaRecoTracksGood", "fLambdaRecoTracksGood", 20, -0.5, 19.5); fLambdaRecoVertices = new TH1F("fLambdaRecoVertices", "fLambdaRecoVertices", 20, -0.5, 19.5); fLambdaMCVertexProton = new TH1F("fLambdaMCVertexProton", "fLambdaMCVertexProton", 20, 0.0, 20.); fLambdaMCVertexAntiProton = new TH1F("fLambdaMCVertexAntiProton", "fLambdaMCVertexAntiProton", 20, 0.0, 20.); fLambdaRecoVertexProton = new TH1F("fLambdaRecoVertexProton", "fLambdaRecoVertexProton", 20, 0.0, 20.); fLambdaRecoVertexAntiProton = new TH1F("fLambdaRecoVertexAntiProton", "fLambdaRecoVertexAntiProton", 20, 0.0, 20.); fLambdaRecoVertexSum = new TH1F("fLambdaRecoVertexSum", "fLambdaRecoVertexSum", 20, 0.0, 20.); fLambdaHitsPerTrackAntiproton = new TH1F("fLambdaHitsPerTrackAntiproton", "fLambdaHitsPerTrackAntiproton", 100, -0.5, 99.5); fLambdaHitsPerTrackPiminus = new TH1F("fLambdaHitsPerTrackPiminus", "fLambdaHitsPerTrackPiminus", 100, -0.5, 99.5); fLambdaHitsPerTrackEplus = new TH1F("fLambdaHitsPerTrackEplus", "fLambdaHitsPerTrackEplus", 100, -0.5, 99.5); fLambdaHitsPerTrackMuplus = new TH1F("fLambdaHitsPerTrackMuplus", "fLambdaHitsPerTrackMuplus", 100, -0.5, 99.5); fLambdaHitsPerTrackProton = new TH1F("fLambdaHitsPerTrackProton", "fLambdaHitsPerTrackProton", 100, -0.5, 99.5); fLambdaHitsPerTrackPiplus = new TH1F("fLambdaHitsPerTrackPiplus", "fLambdaHitsPerTrackPiplus", 100, -0.5, 99.5); fLambdaHitsPerTrackEminus = new TH1F("fLambdaHitsPerTrackEminus", "fLambdaHitsPerTrackEminus", 100, -0.5, 99.5); fLambdaHitsPerTrackMuminus = new TH1F("fLambdaHitsPerTrackMuminus", "fLambdaHitsPerTrackMuminus", 100, -0.5, 99.5); fLambdaHitsPerTrackUnknown = new TH1F("fLambdaHitsPerTrackUnknown", "fLambdaHitsPerTrackUnknown", 100, -0.5, 99.5); fLambdaParticleOriginsAntiproton = new TH1F("fLambdaParticleOriginsAntiproton", "fLambdaParticleOriginsAntiproton", 8000, -4000.5, 3999.5); fLambdaParticleOriginsPiminus = new TH1F("fLambdaParticleOriginsPiminus", "fLambdaParticleOriginsPiminus", 8000, -4000.5, 3999.5); fLambdaParticleOriginsEplus = new TH1F("fLambdaParticleOriginsEplus", "fLambdaParticleOriginsEplus", 8000, -4000.5, 3999.5); fLambdaParticleOriginsMuplus = new TH1F("fLambdaParticleOriginsMuplus", "fLambdaParticleOriginsMuplus", 8000, -4000.5, 3999.5); fLambdaParticleOriginsProton = new TH1F("fLambdaParticleOriginsProton", "fLambdaParticleOriginsProton", 8000, -4000.5, 3999.5); fLambdaParticleOriginsPiplus = new TH1F("fLambdaParticleOriginsPiplus", "fLambdaParticleOriginsPiplus", 8000, -4000.5, 3999.5); fLambdaParticleOriginsEminus = new TH1F("fLambdaParticleOriginsEminus", "fLambdaParticleOriginsEminus", 8000, -4000.5, 3999.5); fLambdaParticleOriginsMuminus = new TH1F("fLambdaParticleOriginsMuminus", "fLambdaParticleOriginsMuminus", 8000, -4000.5, 3999.5); fLambdaParticleOriginsUnknown = new TH1F("fLambdaParticleOriginsUnknown", "fLambdaParticleOriginsUnknown", 8000, -4000.5, 3999.5); fLambdaParticleTypes = new TH1F("fLambdaParticleTypes", "fLambdaParticleTypes", 6000, -3000.5, 2999.5); fLambdaParticleOrigins = new TH1F("fLambdaParticleOrigins", "fLambdaParticleOrigins", 8000, -4000.5, 3999.5); fAlphaResolution = new TH1F("fAlphaResolution", "fAlphaResolution", 1000, -100.00, 100.00); fD1Resolution = new TH1F("fD1Resolution", "fD1Resolution", 1000, -100., 100.); fD0Resolution = new TH1F("fD0Resolution", "fD0Resolution", 1000, -10., 10.); fPhiResolution = new TH1F("fPhiResolution", "fPhiResolution", 1000, -3.2, 3.2); fRadiusResolution = new TH1F("fRadiusResolution", "fRadiusResolution", 1000, -100., 100.); fAlphaReconVsMc = new TH2F("fAlphaResolution", "fAlphaResolution", 100, -200., 200., 1000, -200., 200.); fD1ReconVsMc = new TH2F("fD1Resolution", "fD1Resolution", 100, -400., 400., 400, -400., 400.); fD0ReconVsMc = new TH2F("fD0Resolution", "fD0Resolution", 100, -10., 10., 100, -10., 10.); fPhiReconVsMc = new TH2F("fPhiResolution", "fPhiResolution", 100, -3.2, 3.2, 100, -3.2, 3.2); fRadiusReconVsMc = new TH2F("fRadiusResolution", "fRadiusResolution", 100, -100., 200., 100, -100., 200.); fChiSquareLong = new TH1F("fChiSquareLong", "fChiSquareLong", 1000, -1., 50.); fChiSquareRad = new TH1F("fChiSquareRad","fChiSquareRad", 1000, -1., 50.); fVertexResolutionX = new TH1F("fVertexResolutionX", "fVertexResolutionX", 100, -10., 10.); fVertexResolutionY = new TH1F("fVertexResolutionY", "fVertexResolutionY", 100, -10., 10.); fVertexResolutionZ = new TH1F("fVertexResolutionZ", "fVertexResolutionZ", 100, -10., 10.); fMomentumResolutionX = new TH1F("fMomentumResolutionX", "fMomentumResolutionX", 100, -0.4, 0.4); fMomentumResolutionY = new TH1F("fMomentumResolutionY", "fMomentumResolutionY", 100, -0.4, 0.4); fMomentumResolutionZ = new TH1F("fMomentumResolutionZ", "fMomentumResolutionZ", 100, -0.4, 0.4); fPhiTrResolution = new TH1F("fPhiTrResolution", "fPhiTrResolution", 100, -0.30, 0.30); fThetaResolution = new TH1F("fThetaResolution", "fThetaResolution", 100, -0.30, 0.30); fThetaResolutionCut = new TH1F("fThetaResolutionCut", "fThetaResolutionCut", 100, -0.30, 0.30); fMomentumResolutionTrans = new TH1F("fMomentumResolutionTrans", "fMomentumResolutionTrans", 100, -0.4, 0.4); fMomentumResolutionTransCut = new TH1F("fMomentumResolutionTransCut", "fMomentumResolutionTransCut", 100, -0.4, 0.4); fMomentumResolutionTransPerc = new TH1F("fMomentumResolutionTransPerc", "fMomentumResolutionTransPerc", 100, -2.0, 2.0); fMomentumResolution = new TH1F("fMomentumResolution", "fMomentumResolution", 100, -0.4, 0.4); fVertexDistributionX = new TH2F("fVertexDistributionX", "fVertexDistributionX", 100, -75, 75, 100, -75, 75); fVertexDistributionY = new TH2F("fVertexDistributionY", "fVertexDistributionY", 100, -75, 75, 100, -75, 75); fVertexDistributionZ = new TH2F("fVertexDistributionZ", "fVertexDistributionZ", 100, -75, 75, 100, -75, 75); fMomentumDistributionX = new TH2F("fMomentumDistributionX", "fMomentumDistributionX", 100, -0.1, 1.0, 100, -0.1, 1.0); fMomentumDistributionY = new TH2F("fMomentumDistributionY", "fMomentumDistributionY", 100, -0.1, 1.0, 100, -0.1, 1.0); fMomentumDistributionZ = new TH2F("fMomentumDistributionZ", "fMomentumDistributionZ", 100, -0.1, 1.0, 100, -0.1, 1.0); fPhiTrDistribution = new TH2F("fPhiTrDistribution", "fPhiTrDistribution", 150, -3.20 / 2., 3.20 / 2., 150, -3.20 / 2., 3.20 / 2.); fThetaDistribution = new TH2F("fThetaDistribution", "fThetaDistribution", 150, 0., 3.20 / 2., 150, 0., 3.20 / 2.); fMomentumDistributionTrans = new TH2F("fMomentumDistributionTrans", "fMomentumDistributionTrans", 100, -0.1, 1.0, 100, -0.1, 1.0); fMomentumDistribution = new TH2F("fMomentumDistribution", "fMomentumDistribution", 200, -0.1, 1.5, 200, -0.1, 1.5); fInvariantMass1 = new TH1F("fInvariantMassantiL", "fInvariantMassantiL", 1000, 0., 2.); fInvariantMass2 = new TH1F("fInvariantMassL", "fInvariantMassL", 1000, 0., 2.); fInvMassVsEgy1 = new TH2F("fInvMassVsEgy1", "fInvMassVsEgy1", 1000, 0, 2, 1000, 0, 2); fInvMassVsEgy2 = new TH2F("fInvMassVsEgy2", "fInvMassVsEgy2", 1000, 0, 2, 1000, 0, 2); fE1 = new TH1F("fE1", "fE1", 1000, 0., 10.); fE2 = new TH1F("fE2", "fE2", 1000, 0., 10.); fpx1 = new TH1F("fpx1", "fpx1", 1000, -10., 10.); fpx2 = new TH1F("fpx2", "fpx2", 1000, -10., 10.); fpy1 = new TH1F("fpy1", "fpy1", 1000, -10., 10.); fpy2 = new TH1F("fpy2", "fpy2", 1000, -10., 10.); fpz1 = new TH1F("fpz1", "fpz1", 1000, -5., 5.); fpz2 = new TH1F("fpz2", "fpz2", 1000, -5., 5.); fEnergyAllChiSq = new TH1F("fEnergyAllChiSq", "fEnergyAllChiSq", 1000, 0., 2.); // fParticleIDAllChiSq = new TH1F("fEnergyAllChiSq", "fEnergyAllChiSq", 4600, -2300., 2300.); fParticleIDBadChiSq = new TH1F("fParticleIDBadChiSq", "fParticleIDBadChiSq", 4600, -2300., 2300.); fParticleIDGoodChiSq = new TH1F("fParticleIDGoodChiSq", "fParticleIDGoodChiSq", 4600, -2300., 2300.); fEnergyBadChiSq = new TH1F("fEnergyBadChiSq", "fEnergyBadChiSq", 1000, 0., 2.); fEnergyGoodChiSq = new TH1F("fEnergyGoodChiSq", "fEnergyGoodChiSq", 1000, 0., 2.); fParticleIDGoodChiSqLessThanTF = new TH1F("fParticleIDGoodChiSqLessThanTF", "fParticleIDGoodChiSqLessThanTF", 4600, -2300., 2300. ); fParticleIDGoodChiSqTwentyFiveToThirty = new TH1F("fParticleIDGoodChiSqTwentyFiveToThirty", "fParticleIDGoodChiSqTwentyFiveToThirty", 4600, -2300., 2300.); fParticleIDGoodChiSqSeventyToNinety = new TH1F("fParticleIDGoodChiSqSeventyToNinety", "fParticleIDGoodChiSqSeventyToNinety", 4600, -2300., 2300.); fParticleIDBadChiSqLessThanTh = new TH1F("fParticleIDBadChiSqLessThanTh", "fParticleIDBadChiSqLessThanTh", 4600, -2300., 2300. ); fParticleIDBadChiSqNineToEleven = new TH1F("fParticleIDBadChiSqNineToEleven", "fParticleIDBadChiSqNineToEleven", 4600, -2300., 2300.); fParticleIDBadChiSqFifteenToNineteen = new TH1F("fParticleIDBadChiSqFifteenToNineteen", "fParticleIDBadChiSqFifteenToNineteen", 4600, -2300., 2300.); fXVertexMcMinReco = new TH1F("fXVertexMcMinReco","fXVertexMcMinReco", 1000, -10., 10.); fYVertexMcMinReco = new TH1F("fYVertexMcMinReco","fYVertexMcMinReco", 1000, -10., 10.); fZVertexMcMinReco = new TH1F("fZVertexMcMinReco","fZVertexMcMinReco", 1000, -10., 10.); } void CbmSttFitTracksQa::WriteHistograms(char *filename) { TFile *hfile = new TFile(filename,"RECREATE","histograms for CbmSttFitTracksQa"); fLambdaVertexMomentumAngleLambda->Write(); fLambdaVertexMomentumAngleAntiLambda->Write(); fLambdaBadPBarTheta->Write(); fLambdaGoodPBarTheta->Write(); fLambdaBadPBarHits->Write(); fLambdaGoodPBarHits->Write(); fLambdaSelectorHits->Write(); fLambdaSelectorHitsBad->Write(); fLambdaProtonChiSquareLong->Write(); fLambdaPiPlusChiSquareLong->Write(); fLambdaPiMinusChiSquareLong->Write(); fLambdaAntiProtonChiSquareLong->Write(); fLambdaProtonChiSquareRad->Write(); fLambdaPiPlusChiSquareRad->Write(); fLambdaPiMinusChiSquareRad->Write(); fLambdaAntiProtonChiSquareRad->Write(); fLambdaCosThetaLambdaMc->Write(); fLambdaCosThetaLambdaReco->Write(); fLambdaCosThetaAntiLambdaReco->Write(); fLambdaPhiLambdaReco->Write(); fLambdaPhiAntiLambdaReco->Write(); fLambdaOpeningPhiLambdaLambdabarReco->Write(); fLambdaOpeningThetaLambdaLambdabarReco->Write(); fLambdaTotalPxReco->Write(); fLambdaTotalPyReco->Write(); fLambdaTotalPzReco->Write(); fLambdaTotalEnReco->Write(); fLambdaTotalPxMc->Write(); fLambdaTotalPyMc->Write(); fLambdaTotalPzMc->Write(); fLambdaTotalEnMc->Write(); fLambdaInvMassLambda->Write(); fLambdaInvMassAntiLambda->Write(); fLambdaInvMassSum->Write(); fLambdaLambdaMomentumX->Write(); fLambdaLambdaMomentumY->Write(); fLambdaLambdaMomentumZ->Write(); fLambdaAntiLambdaMomentumX->Write(); fLambdaAntiLambdaMomentumY->Write(); fLambdaAntiLambdaMomentumZ->Write(); fLambdaPiPlusMomentumX->Write(); fLambdaPiPlusMomentumY->Write(); fLambdaPiPlusMomentumZ->Write(); fLambdaPiPlusMomentumPt->Write(); fLambdaPiPlusMomentumP->Write(); fLambdaPiMinusMomentumX->Write(); fLambdaPiMinusMomentumY->Write(); fLambdaPiMinusMomentumZ->Write(); fLambdaPiMinusMomentumPt->Write(); fLambdaPiMinusMomentumP->Write(); fLambdaProtonMomentumX->Write(); fLambdaProtonMomentumY->Write(); fLambdaProtonMomentumZ->Write(); fLambdaProtonMomentumPt->Write(); fLambdaProtonMomentumP->Write(); fLambdaAntiProtonMomentumX->Write(); fLambdaAntiProtonMomentumY->Write(); fLambdaAntiProtonMomentumZ->Write(); fLambdaAntiProtonMomentumPt->Write(); fLambdaAntiProtonMomentumP->Write(); fLambdaVertexResolutionX->Write(); fLambdaVertexResolutionY->Write(); fLambdaVertexResolutionZ->Write(); fLambdaVertexResolution->Write(); fLambdaRecoVertices->Write(); fLambdaMCVertexProton->Write(); fLambdaMCVertexAntiProton->Write(); fLambdaRecoVertexProton->Write(); fLambdaRecoVertexAntiProton->Write(); fLambdaRecoVertexSum->Write(); fLambdaRecoTracks->Write(); fLambdaParticleTypes->Write(); fLambdaParticleOrigins->Write(); fLambdaRecoTracksAbove->Write(); fLambdaRecoTracksAnnihilation->Write(); fLambdaRecoTracksPionDecay->Write(); fLambdaRecoTracksGood->Write(); fLambdaHitsPerTrack->Write(); fLambdaHitsPerTrackAntiproton->Write(); fLambdaHitsPerTrackPiminus->Write(); fLambdaHitsPerTrackEplus->Write(); fLambdaHitsPerTrackMuplus->Write(); fLambdaHitsPerTrackProton->Write(); fLambdaHitsPerTrackPiplus->Write(); fLambdaHitsPerTrackEminus->Write(); fLambdaHitsPerTrackMuminus->Write(); fLambdaHitsPerTrackUnknown->Write(); fLambdaParticleOriginsAntiproton->Write(); fLambdaParticleOriginsPiminus->Write(); fLambdaParticleOriginsEplus->Write(); fLambdaParticleOriginsMuplus->Write(); fLambdaParticleOriginsProton->Write(); fLambdaParticleOriginsPiplus->Write(); fLambdaParticleOriginsEminus->Write(); fLambdaParticleOriginsMuminus->Write(); fLambdaParticleOriginsUnknown->Write(); fAlphaResolution->Write(); fD1Resolution->Write(); fD0Resolution->Write(); fPhiResolution->Write(); fRadiusResolution->Write(); fAlphaReconVsMc->Write(); fD1ReconVsMc->Write(); fD0ReconVsMc->Write(); fPhiReconVsMc->Write(); fRadiusReconVsMc->Write(); fChiSquareLong->Write(); fChiSquareRad->Write(); fVertexResolutionX->Write(); fVertexResolutionY->Write(); fVertexResolutionZ->Write(); fMomentumResolutionX->Write(); fMomentumResolutionY->Write(); fMomentumResolutionZ->Write(); fPhiTrResolution->Write(); fThetaResolution->Write(); fThetaResolutionCut->Write(); fMomentumResolutionTrans->Write(); fMomentumResolutionTransCut->Write(); fMomentumResolutionTransPerc->Write(); fMomentumResolution->Write(); fVertexDistributionX->Write(); fVertexDistributionY->Write(); fVertexDistributionZ->Write(); fMomentumDistributionX->Write(); fMomentumDistributionY->Write(); fMomentumDistributionZ->Write(); fPhiTrDistribution->Write(); fThetaDistribution->Write(); fMomentumDistributionTrans->Write(); fMomentumDistribution->Write(); fInvariantMass1->Write(); fInvariantMass2->Write(); fInvMassVsEgy1->Write(); fInvMassVsEgy2->Write(); fE1->Write(); fE2->Write(); fpx1->Write(); fpx2->Write(); fpy1->Write(); fpy2->Write(); fpz1->Write(); fpz2->Write(); fParticleIDBadChiSq->Write(); fParticleIDGoodChiSq->Write(); fEnergyBadChiSq->Write(); fEnergyGoodChiSq->Write(); fParticleIDGoodChiSqLessThanTF->Write(); fParticleIDGoodChiSqTwentyFiveToThirty->Write(); fParticleIDGoodChiSqSeventyToNinety->Write(); fEnergyAllChiSq->Write(); fParticleIDBadChiSqLessThanTh->Write(); fParticleIDBadChiSqNineToEleven->Write(); fParticleIDBadChiSqFifteenToNineteen->Write(); fXVertexMcMinReco->Write(); fYVertexMcMinReco->Write(); fZVertexMcMinReco->Write(); hfile->Close(); delete hfile; } void CbmSttFitTracksQa::drawDetector() { Double_t InnerRadius = 16., OuterRadius = 42.; Float_t pointsxo[100], pointsyo[100], pointsxi[100], pointsyi[100], pointszb[100], pointsza[100]; for (Int_t teller = 0; teller < 100; teller++) { pointsxo[teller] = OuterRadius * cos(teller * (2 * dPi / 100)); pointsyo[teller] = OuterRadius * sin(teller * (2 * dPi / 100)); pointsxi[teller] = InnerRadius * cos(teller * (2 * dPi / 100)); pointsyi[teller] = InnerRadius * sin(teller * (2 * dPi / 100)); pointszb[teller] = -75.; pointsza[teller] = 75.; } TPolyLine3D *myLineAO = new TPolyLine3D(100, pointsxo, pointsyo, pointsza); TPolyLine3D *myLineBO = new TPolyLine3D(100, pointsxo, pointsyo, pointszb); TPolyLine3D *myLineAI = new TPolyLine3D(100, pointsxi, pointsyi, pointsza); TPolyLine3D *myLineBI = new TPolyLine3D(100, pointsxi, pointsyi, pointszb); myLineAO->Draw(); myLineBO->Draw(); myLineAI->Draw(); myLineBI->Draw(); } void CbmSttFitTracksQa::drawAxis() { Float_t xaxisx[2] = {-rangex, rangex}, yaxisx[2] = {0.0, 0.0}, zaxisx[2] = {0.0, 0.0}; Float_t xaxisy[2] = {0.0, 0.0}, yaxisy[2] = {-rangey, rangey}, zaxisy[2] = {0.0, 0.0}; Float_t xaxisz[2] = {0.0, 0.0}, yaxisz[2] = {0.0, 0.0}, zaxisz[2] = {-rangez, rangez}; TPolyLine3D *xaxis = new TPolyLine3D(2, xaxisx, xaxisy, xaxisz), *yaxis = new TPolyLine3D(2, yaxisx, yaxisy, yaxisz), *zaxis = new TPolyLine3D(2, zaxisx, zaxisy, zaxisz); xaxis->SetLineColor(1); xaxis->Draw(); yaxis->SetLineColor(2); yaxis->Draw(); zaxis->SetLineColor(3); zaxis->Draw(); } void CbmSttFitTracksQa::drawHits(Int_t trackno) { CbmSttTrack *track = (CbmSttTrack *) fSttTracks->At(trackno); for (Int_t hitCounter = 0; hitCounter < track->GetNofHits(); hitCounter++) { CbmSttHit *hit = (CbmSttHit *) fSttHits->At(track->GetHitIndex(hitCounter)); CbmSttGeomPoint myPoint(hit->GetX(), hit->GetY(), hit->GetZ()); myPoint.Draw(1.0); } } Int_t CbmSttFitTracksQa::LambdaLambdabarSelector(CbmSttTrack *&proton, CbmSttTrack *&negpion, CbmSttTrack *&antiproton, CbmSttTrack *&pospion) { CbmSttVertex *lambdaVertex, *antilambdaVertex; Int_t retval = 0; Int_t topoTeller = 0; CbmSttTrack *bestProton, *bestAntiProton, *bestPosPion, *bestNegPion; Double_t bestEnergyBalance = 10000.; while (GetTopology(lambdaVertex, antilambdaVertex, proton, negpion, antiproton, pospion, topoTeller)) { Double_t energyBalance; if (CheckKinematics(lambdaVertex, antilambdaVertex, proton, negpion, antiproton, pospion, energyBalance)) { retval++; if (energyBalance < bestEnergyBalance) { bestProton = proton; bestAntiProton = antiproton; bestPosPion = pospion; bestNegPion = negpion; bestEnergyBalance = energyBalance; } } topoTeller++; } proton = bestProton; antiproton = bestAntiProton; pospion = bestPosPion; negpion = bestNegPion; return retval; } void CbmSttFitTracksQa::ConnectTracksToVertices(vector &myTracksVertex1, vector &myTracksVertex2) { // if 2 vertices if (fSttVertices->GetEntriesFast() == 2) { Int_t nRec = fSttTracks->GetEntriesFast(); // loop over all reconstructed particles for (Int_t iRec=0; iRec < nRec; iRec++) { CbmSttTrack *track = (CbmSttTrack*) fSttTracks->At(iRec); if (!track) { cout << "-W- " << GetName() << "::Exec: Empty SttTrack at " << iRec << endl; continue; } // if > 4 hits and charge determined, put in array if (track->GetNHits() > 4) { if (fabs(track->GetParamLast()->GetQp()) > 0.) { Double_t tmp1, tmp2, tmp3; if (GetVertex(track, tmp1, tmp2, tmp3) == 0) { //cout << "adding to 1" << endl; myTracksVertex1.push_back(track); } else { if (GetVertex(track, tmp1, tmp2, tmp3) == 1) { //cout << "adding to 2" << endl; myTracksVertex2.push_back(track); } } } } } } } Bool_t CbmSttFitTracksQa::GetTopology(CbmSttVertex *&lambdaVertex, CbmSttVertex *&antilambdaVertex, CbmSttTrack *&proton, CbmSttTrack *&negpion, CbmSttTrack *&antiproton, CbmSttTrack *&pospion, Int_t index) { Bool_t retval = kFALSE; Int_t nUsable = 0, nRec = fSttTracks->GetEntriesFast(); // loop over all reconstructed particles for (Int_t iRec=0; iRec < nRec; iRec++) { CbmSttTrack *track = (CbmSttTrack*) fSttTracks->At(iRec); if (!track) { cout << "-W- " << GetName() << "::Exec: Empty SttTrack at " << iRec << endl; continue; } if (track->GetNHits() > 4) { nUsable++; } } if (nUsable == 4) { // TODO: check chisquare on vertices // if 2 vertices if (fSttVertices->GetEntriesFast() == 2) { vector myTracksVertex1, myTracksVertex2; ConnectTracksToVertices(myTracksVertex1, myTracksVertex2); // if 2 particles per vertex if ((myTracksVertex1.size() == 2) && (myTracksVertex2.size() == 2)) { // positive track and negative track at each vertex if ((myTracksVertex1[0]->GetParamLast()->GetQp() * myTracksVertex1[1]->GetParamLast()->GetQp() < 0.) && (myTracksVertex2[0]->GetParamLast()->GetQp() * myTracksVertex2[1]->GetParamLast()->GetQp() < 0.)) { // two possible topologies lambda = vertex1 and antilambda = vertex2 or vice versa if (index == 0) { if (myTracksVertex1[0]->GetParamLast()->GetQp() > 0.) { proton = myTracksVertex1[0]; negpion = myTracksVertex1[1]; } else { proton = myTracksVertex1[1]; negpion = myTracksVertex1[0]; } if (myTracksVertex2[0]->GetParamLast()->GetQp() > 0.) { pospion = myTracksVertex2[0]; antiproton = myTracksVertex2[1]; } else { pospion = myTracksVertex2[1]; antiproton = myTracksVertex2[0]; } antilambdaVertex = (CbmSttVertex*) fSttVertices->At(1); lambdaVertex = (CbmSttVertex*) fSttVertices->At(0); retval = kTRUE; } else if (index == 1) { if (myTracksVertex1[0]->GetParamLast()->GetQp() > 0.) { antiproton = myTracksVertex1[1]; pospion = myTracksVertex1[0]; } else { antiproton = myTracksVertex1[0]; pospion = myTracksVertex1[1]; } if (myTracksVertex2[0]->GetParamLast()->GetQp() > 0.) { proton = myTracksVertex2[0]; negpion = myTracksVertex2[1]; } else { proton = myTracksVertex2[1]; negpion = myTracksVertex2[0]; } antilambdaVertex = (CbmSttVertex*) fSttVertices->At(0); lambdaVertex = (CbmSttVertex*) fSttVertices->At(1); retval = kTRUE; } } } } } return retval; } Bool_t CbmSttFitTracksQa::CheckKinematics(CbmSttVertex *lambdaVertex, CbmSttVertex *antilambdaVertex, CbmSttTrack *proton, CbmSttTrack *negpion, CbmSttTrack *antiproton, CbmSttTrack *pospion, Double_t &energyBal) { Bool_t retval = kFALSE; // get helices CbmSttGeomHelix protonHel(proton), negpionHel(negpion), antiprotonHel(antiproton), pospionHel(pospion); // get vertex CbmSttGeomPoint lambdaVertexPoint(lambdaVertex->GetVertexX(), lambdaVertex->GetVertexY(), lambdaVertex->GetVertexZ()), antilambdaVertexPoint(antilambdaVertex->GetVertexX(), antilambdaVertex->GetVertexY(), antilambdaVertex->GetVertexZ()); TVector3 vertexvecLambda(lambdaVertex->GetVertexX(), lambdaVertex->GetVertexY(), lambdaVertex->GetVertexZ()), vertexvecAntiLambda(antilambdaVertex->GetVertexX(), antilambdaVertex->GetVertexY(), antilambdaVertex->GetVertexZ()); // get masses Double_t massProton = TDatabasePDG::Instance()->GetParticle(2212)->Mass(), massAntiProton = TDatabasePDG::Instance()->GetParticle(-2212)->Mass(), massPiPlus = TDatabasePDG::Instance()->GetParticle(211)->Mass(), massPiMinus = TDatabasePDG::Instance()->GetParticle(-211)->Mass(); // get energy Double_t protonEn = sqrt(protonHel.GetP() * protonHel.GetP() + massProton * massProton), negpionEn = sqrt(negpionHel.GetP() * negpionHel.GetP() + massPiMinus * massPiMinus), antiprotonEn = sqrt(antiprotonHel.GetP() * antiprotonHel.GetP() + massAntiProton * massAntiProton), pospionEn = sqrt(pospionHel.GetP() * pospionHel.GetP() + massPiPlus * massPiPlus); // get momentum TLorentzVector protonVec(protonHel.GetMomentumAt(lambdaVertexPoint), protonEn), negpionVec(negpionHel.GetMomentumAt(lambdaVertexPoint), negpionEn), antiprotonVec(antiprotonHel.GetMomentumAt(antilambdaVertexPoint), antiprotonEn), pospionVec(pospionHel.GetMomentumAt(antilambdaVertexPoint), pospionEn); TLorentzVector lambdaVec = protonVec + negpionVec, antilambdaVec = antiprotonVec + pospionVec; Double_t lambdaAngle = vertexvecLambda.Angle(lambdaVec.Vect()), antilambdaAngle = vertexvecAntiLambda.Angle(antilambdaVec.Vect()); TLorentzVector initialState(0., 0., 2., massProton + sqrt(massAntiProton * massAntiProton + 2. * 2.)), total = lambdaVec + antilambdaVec - initialState; // total energy balance at 0 if (fabs(total.Energy()) < 3 * 0.17) { retval = kTRUE; } energyBal = fabs(total.Energy()); return retval; } Int_t CbmSttFitTracksQa::GetVertex(CbmSttTrack *recoTrack, Double_t &recoVertexProtonX, Double_t &recoVertexProtonY, Double_t &recoVertexProtonZ) { Int_t retval = -1; //cout << "looking for track: " << recoTrack << endl; // loop over all vertices Int_t nVertices = fSttVertices->GetEntriesFast(); for (Int_t vertCounter = 0; vertCounter < nVertices; vertCounter++) { CbmSttVertex *vertex = (CbmSttVertex*) fSttVertices->At(vertCounter); //cout << "vertex: " << vertCounter << " track1: " << (CbmSttTrack *)fSttTracks->At(vertex->GetTrack(0)) << endl; //cout << "vertex: " << vertCounter << " track2: " << (CbmSttTrack *)fSttTracks->At(vertex->GetTrack(1)) << endl; if (((CbmSttTrack *)fSttTracks->At(vertex->GetTrack(0)) == recoTrack) || ((CbmSttTrack *)fSttTracks->At(vertex->GetTrack(1)) == recoTrack)) { retval = vertCounter; recoVertexProtonX = vertex->GetVertexX(); recoVertexProtonY = vertex->GetVertexY(); recoVertexProtonZ = vertex->GetVertexZ(); } } return retval; } ClassImp(CbmSttFitTracksQa)