/* $Id: CbmRichAnalysisHits.cxx,v 1.7 2006/02/09 16:30:46 hoehne Exp $ */ /* History of cvs commits: * * $Log: CbmRichAnalysisHits.cxx,v $ * Revision 1.7 2006/02/09 16:30:46 hoehne * changes due to removed trackID in CbmRichHit * * Revision 1.6 2006/01/11 12:42:40 hoehne * some usage of CbmRichPhotodetector replaced * * Revision 1.5 2005/12/19 19:04:30 friese * New CbmTask design * * Revision 1.4 2005/09/08 15:43:29 kharlov * Set flag to analysis types * * Revision 1.3 2005/07/18 09:55:08 kharlov * Track projection analysis is added * * Revision 1.2 2005/07/06 14:59:45 kharlov * Non-used histos are removed * */ #include "CbmRichAnalysisHits.h" #include "CbmRootManager.h" #include "CbmMCTrack.h" #include "CbmRichPoint.h" #include "CbmRichHit.h" #include "CbmRichRing.h" #include "CbmRichLightSpot.h" #include "CbmRichLightSpotMC.h" #include "CbmRichProjectionProducer.h" #include "iostream.h" #include "TParticle.h" #include "TH1F.h" #include "TH2F.h" #include "TH3F.h" #include "TGraph.h" #include "TVector3.h" //------------------------------------------------------------------------------ CbmRichAnalysisHits::CbmRichAnalysisHits() :CbmTask() { } //------------------------------------------------------------------------------ CbmRichAnalysisHits::CbmRichAnalysisHits(const char *name, const char *title) :CbmTask(name) { fEvent = 0; SetDebug(); SetWhatToAnalyze(); SetEbyEgraph(); const Int_t maxHit = 5000; const Int_t maxRG = 1000; fNHits = 0; fXHits = new TArrayD(maxHit); fYHits = new TArrayD(maxHit); fNrgs = 0; fXrgs = new TArrayD(maxRG); fYrgs = new TArrayD(maxRG); fXtr = new TArrayD(maxRG); fYtr = new TArrayD(maxRG); fXcp = new TArrayD(maxRG); fYcp = new TArrayD(maxRG); Double_t pi = TMath::Pi(); // Primary tracks histograms if (strstr(fWhatToAnalyze,"T")) { fhPrimPabs = new TH1F("hPrimPabs" ,"Momentum of primary particles",100, 0.,50.); fhPrimPt = new TH1F("hPrimPt" ,"p_{T} of primary particles" ,100, 0.,10.); fhPrimEta = new TH1F("hPrimEta" ,"#eta of primary particles" , 80,-2., 6.); fhPrimTheta= new TH1F("hPrimTheta","#theta of primary particles" ,100, 0., pi); fhPrimPhi = new TH1F("hPrimPhi" ,"#phi of primary particles" ,100, 0.,2*pi); fhPrimPtEta= new TH2F("hPrimPtEta","p_{T},#eta of primary particles",100, 0.,10., 80,-2., 6.); } // MC points histograms if (strstr(fWhatToAnalyze,"P")) { fhNPhotonsPrim = new TH1F("hNPhotonsPrim","Number of fired PMTs",300,0.,300.); fhPhotonsXYmcp = new TH2F("hPhotonsXYmcp","(X,Y) coordinates of Cherenkov photons", 500,-163,163,850,-255,255); fhPhotonsE = new TH1F("hPhotonsE","Cherenkov photons energy",100,0.,10.); fhEtaVsX = new TH2F("hEtaVsX","#eta vs Cherenkov X",40,-2.,6.,50,0.,163); fhEtaVsY = new TH2F("hEtaVsY","#eta vs Cherenkov Y",40,-2.,6.,50,0.,255); fhEtaVsR = new TH2F("hEtaVsR","#eta vs Cherenkov R",40,-2.,6.,50,0.,300); } // Hits histograms if (strstr(fWhatToAnalyze,"H")) { fhNpmt = new TH1F("hNpmt" ,"Number of fires PMTs",500,0.,maxHit); fhPhotonsXY = new TH2F("hPhotonsXY","(X,Y) coordinates of Cherenkov photons", 500,-163,163,850,-255,255); fhNTracks = new TH1F("hNTracks" ,"Number of detected tracks",50,0.,200.); fhNPrimTracks = new TH1F("hNPrimTracks" ,"Number of detected primary tracks",50,0.,50.); fhTrackZ = new TH1F("hTrackZ" ,"Z of detected tracks", 400,0.,400.); fhNphotPMT = new TH1F("hNphotPMT","Number of photons per PMT",10,0.5,10.5); fhAmplitude = new TH1F("hAmplitude","Signal amplitude in PMT",100,0.,10.); } // Kinematics of detected tracks per their type if (strstr(fWhatToAnalyze,"R")) { fhPtElectron = new TH1F("hPtElectron" ,"p_{T} of detected e^{#pm}" ,50,0.,5.); fhPtMuon = new TH1F("hPtMuon" ,"p_{T} of detected #mu^{#pm}",50,0.,5.); fhPtPion = new TH1F("hPtPion" ,"p_{T} of detected #pi^{#pm}",50,0.,5.); fhPtKaon = new TH1F("hPtKaon" ,"p_{T} of detected K^{#pm}" ,50,0.,5.); fhPtProton = new TH1F("hPtProton" ,"p_{T} of detected p^{#pm}" ,50,0.,5.); fhEtaElectron = new TH1F("hEtaElectron","#eta of detected e^{#pm}" ,50,0.,5.); fhEtaMuon = new TH1F("hEtaMuon" ,"#eta of detected #mu^{#pm}" ,50,0.,5.); fhEtaPion = new TH1F("hEtaPion" ,"#eta of detected #pi^{#pm}" ,50,0.,5.); fhEtaKaon = new TH1F("hEtaKaon" ,"#eta of detected K^{#pm}" ,50,0.,5.); fhEtaProton = new TH1F("hEtaProton" ,"#eta of detected p^{#pm}" ,50,0.,5.); fhPElectron = new TH1F("hPElectron" ,"p of detected e^{#pm}" ,300,0.,30.); fhPMuon = new TH1F("hPMuon" ,"p of detected #mu^{#pm}" ,300,0.,30.); fhPPion = new TH1F("hPPion" ,"p of detected #pi^{#pm}" ,300,0.,30.); fhPKaon = new TH1F("hPKaon" ,"p of detected K^{#pm}" ,300,0.,30.); fhPProton = new TH1F("hPProton" ,"p of detected p^{#pm}" ,300,0.,30.); // Light spot characteristics fhNLSElectron = new TH2F("hNLSElectron","N_{PMT} vs P for e^{#pm}" ,75,0.,15.,100,0.5,100.5); fhRLSElectron = new TH2F("hRLSElectron","R_{ring} vs P for e^{#pm}" ,75,0.,15.,100,0.,10.); fhNLSMuon = new TH2F("hNLSMuon" ,"N_{PMT} vs P for #mu^{#pm}" ,75,0.,15.,100,0.5,100.5); fhRLSMuon = new TH2F("hRLSMuon" ,"R_{ring} vs P for #mu^{#pm}",75,0.,15.,100,0.,10.); fhNLSPion = new TH2F("hNLSPion" ,"N_{PMT} vs P for #pi^{#pm}" ,75,0.,15.,100,0.5,100.5); fhRLSPion = new TH2F("hRLSPion" ,"R_{ring} vs P for #pi^{#pm}",75,0.,15.,100,0.,10.); fhNLSKaon = new TH2F("hNLSKaon" ,"N_{PMT} vs P for K^{#pm}" ,75,0.,15.,100,0.5,100.5); fhRLSKaon = new TH2F("hRLSKaon" ,"R_{ring} vs P for K^{#pm}" ,75,0.,15.,100,0.,10.); fhNLSProton = new TH2F("hNLSProton" ,"N_{PMT} vs P for p^{#pm}" ,75,0.,15.,100,0.5,100.5); fhRLSProton = new TH2F("hRLSProton" ,"R_{ring} vs P for p^{#pm}" ,75,0.,15.,100,0.,10.); // Distance between light spot center and ring guidance fhRtrueMatch = new TH1F("hRtrueMatch" ,"#Delta_{true} (ls-rg)" ,500,0.,100.); fhRfalseMatch = new TH1F("hRfalseMatch","#Delta_{false} (ls-rg)",500,0.,100.); // Distance between light spot center and ring guidance for closest distance match fhTrue = new TH1F("hTrue" ,"#Delta_{true} (ls-rg)" ,500,0.,100.); fhTrueElectron = new TH1F("hTrueElectron" ,"#Delta_{true} (ls-rg) for e" ,500,0.,100.); fhTruePion = new TH1F("hTruePion" ,"#Delta_{true} (ls-rg) for #pi" ,500,0.,100.); fhFalse = new TH1F("hFalse" ,"#Delta_{false} (ls-rg)" ,500,0.,100.); fhFalseElectronPi = new TH1F("hFalseElectronPi" ,"#Delta_{false} (ls-rg) for pi id as e" ,500,0.,100.); fhFalseElectron = new TH1F("hFalseElectron" ,"#Delta_{false} (ls-rg) for any id as e" ,500,0.,100.); // Normalized distance between light spot center and ring guidance fhPdfElectron=new TH1F("hPdfElectron","PDF of #Delta_{true}(ls-rg) for e",100,0.,10.); fhPdfPi = new TH1F("hPdfPi","PDF of #Delta_{true}(ls-rg) for #pi",100,0.,10.); } // Acceptance if (strstr(fWhatToAnalyze,"A")) { fhPtYrap = new TH2F("PtYrap","pt vs y of primary particles", 100,0.,5.,100,0.,6.); fhPtYrapSignal = new TH2F("PtYrapSignal","pt vs y of signal particles", 100,0.,5.,100,0.,6.); fhPtYrapAccept = new TH2F("PtYrapAccept","RICH p_{T}-y acceptance",100,0.,5.,100,0.,6.); // Invariant mass of dilepton pair fhDileptonMass = new TH1F("DileptonMass","Invariant mass of registered electrons",100,0.,4.); fhNphotMCRing = new TH1F("fhNphotMCRing","Number of photons per MC ring",300,0.,300.); } } //------------------------------------------------------------------------------ CbmRichAnalysisHits::~CbmRichAnalysisHits() { // Primary tracks histograms if (strstr(fWhatToAnalyze,"T")) { fhPrimPabs ->Write(); fhPrimPt ->Write(); fhPrimEta ->Write(); fhPrimTheta->Write(); fhPrimPhi ->Write(); fhPrimPtEta->Write(); } // MC points histograms if (strstr(fWhatToAnalyze,"P")) { fhNPhotonsPrim->Write(); fhPhotonsXYmcp->Write(); fhPhotonsE ->Write(); fhEtaVsX ->Write(); fhEtaVsY ->Write(); fhEtaVsR ->Write(); } // Hits histograms if (strstr(fWhatToAnalyze,"H")) { fhNpmt ->Write(); fhPhotonsXY ->Write(); fhNTracks ->Write(); fhNPrimTracks ->Write(); fhTrackZ ->Write(); fhNphotPMT ->Write(); fhAmplitude ->Write(); } // Kinematics of detected tracks per their type if (strstr(fWhatToAnalyze,"R")) { fhPtElectron ->Write(); fhPtMuon ->Write(); fhPtPion ->Write(); fhPtKaon ->Write(); fhPtProton ->Write(); fhEtaElectron ->Write(); fhEtaMuon ->Write(); fhEtaPion ->Write(); fhEtaKaon ->Write(); fhEtaProton ->Write(); fhPElectron ->Write(); fhPMuon ->Write(); fhPPion ->Write(); fhPKaon ->Write(); fhPProton ->Write(); fhNLSElectron ->Write(); fhRLSElectron ->Write(); fhNLSMuon ->Write(); fhRLSMuon ->Write(); fhNLSPion ->Write(); fhRLSPion ->Write(); fhNLSKaon ->Write(); fhRLSKaon ->Write(); fhNLSProton ->Write(); fhRLSProton ->Write(); fhRtrueMatch ->Write(); fhRfalseMatch ->Write(); // Distance between light spot center and ring guidance for closest distance match fhTrue ->Write(); fhTrueElectron ->Write(); fhTruePion ->Write(); fhFalse ->Write(); fhFalseElectronPi ->Write(); fhFalseElectron ->Write(); fhPdfElectron->Sumw2(); fhPdfElectron->Scale(1./fhPdfElectron->GetEntries()); fhPdfElectron ->Write(); fhPdfPi->Sumw2(); fhPdfPi->Scale(1./fhPdfPi->GetEntries()); fhPdfPi ->Write(); //Ring radius resolution // fhRadMCXY->Write(); // fhRadHitsXY->Write(); fhNphotMCRing->Write(); } //Acceptance histograms if (strstr(fWhatToAnalyze,"A")) { fhPtYrap->Sumw2(); fhPtYrap->Write(); fhPtYrapSignal->Sumw2(); fhPtYrapSignal->Write(); fhPtYrapAccept->Divide(fhPtYrapSignal,fhPtYrap); fhPtYrapAccept->Write(); fhDileptonMass->Write(); } } //------------------------------------------------------------------------------ InitStatus CbmRichAnalysisHits::Init() { CbmRootManager *manager= CbmRootManager::Instance(); // all tracks fListStack = (TClonesArray*)manager->ActivateBranch("MCTrack"); //RICH MC points if (strstr(fWhatToAnalyze,"P")) fListRICHpts = (TClonesArray*)manager->ActivateBranch("RichPoint"); //RICH hits if (strstr(fWhatToAnalyze,"H")) fListRICHhits = (TClonesArray*)manager->ActivateBranch("RichHit"); if (strstr(fWhatToAnalyze,"R")) { //RICH ring guidances fListRICHrings = (TClonesArray*)manager->ActivateBranch("RichRingGuidance"); //RICH track projections fProjectionCollection = (TClonesArray*)manager->ActivateBranch("RichProjection"); //RICH track extrapolations fTrackCollection = (TClonesArray*)manager->ActivateBranch("RichTrack"); // light spots fListLightSpots = new TClonesArray("CbmRichLightSpot",500); } } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::Exec(Option_t* option) { // make RICH hit analysis printf("RICH hit analysis: event %d\n",fEvent); if (strstr(fWhatToAnalyze,"T")) AnalyzePrimary(); if (strstr(fWhatToAnalyze,"P")) AnalyzeMCPoints(); if (strstr(fWhatToAnalyze,"H")) AnalyzeHits(); if (strstr(fWhatToAnalyze,"R")) { AnalyzeRingGuidances(); AnalyzeRings(); } if (strstr(fWhatToAnalyze,"A")) { AnalyzeAcceptanceElectrons(); AnalyzeAcceptanceVectorMesons(); } fEvent++; } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::AnalyzePrimary() { // Primary track analysis CbmMCTrack *primary; fNprimary = 0; TLorentzVector momentum; for (Int_t iPrimary=0; iPrimaryGetEntries(); iPrimary++) { primary = (CbmMCTrack*)fListStack->At(iPrimary); if (primary->GetMotherID() != -1) continue; fNprimary++; momentum = primary->Get4Momentum(); if (strstr(fDebug,"prim")) printf("track %d: id=%d, p=(%f,%f,%f,%f) GeV/c, pt,y,theta=%f,%f,%f\n", iPrimary, primary->GetPdgCode(), momentum.Px(), momentum.Py(), momentum.Pz(), momentum.Energy(), momentum.Pt(), momentum.Eta(), momentum.Theta()); fhPrimPabs ->Fill(momentum.P()); fhPrimPt ->Fill(momentum.Pt()); fhPrimEta ->Fill(momentum.Eta()); fhPrimTheta->Fill(momentum.Theta()); fhPrimPhi ->Fill(momentum.Phi()); fhPrimPtEta->Fill(momentum.Pt(),momentum.Eta()); } if (strstr(fDebug,"prim")) printf("\tNumber of primaries: %d\n",fNprimary); } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::AnalyzeMCPoints() { // Monte-Carlo points analysis if (strstr(fDebug,"mcp")) printf("\tNumber of RICH MC points: %d\n",fListRICHpts->GetEntries()); CbmRichPoint* mcPoint=NULL; Int_t nPhotonsPrimary = 0; TVector3 xyz, pxyz; for (Int_t iPoint=0; iPointGetEntries(); iPoint++) { mcPoint = (CbmRichPoint*)fListRICHpts->At(iPoint); // MC point CbmMCTrack* cherenkoff = (CbmMCTrack*)fListStack->At(mcPoint->GetTrackID()); if (cherenkoff->GetPdgCode() != 50000050 ) continue; // not a Cherenkov photon! mcPoint->Position(xyz); mcPoint->Momentum(pxyz); pxyz *= 1e+9; Int_t mother = cherenkoff->GetMotherID(); fhPhotonsXYmcp->Fill(xyz.X(),xyz.Y()); fhPhotonsE ->Fill(pxyz.Mag()); // Analyze only the very primary electron if (mother != 0) continue; nPhotonsPrimary++; // YVK! CbmMCTrack *primaryElectron = (CbmMCTrack*)fListStack->At(mother); Float_t etaMother = primaryElectron->Get4Momentum().Eta(); fhEtaVsX ->Fill(etaMother,TMath::Abs(xyz.X())); fhEtaVsY ->Fill(etaMother,TMath::Abs(xyz.Y())); fhEtaVsR ->Fill(etaMother,xyz.Perp()); } if (nPhotonsPrimary > 0) fhNPhotonsPrim->Fill(nPhotonsPrimary); } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::AnalyzeHits() { // Hit analysis if (strstr(fDebug,"hit")) printf("\tNumber of RICH hits: %d\n",fListRICHhits->GetEntries()); CbmRichHit* hit=NULL; Int_t nLightSpots = 0; fNHits = 0; Int_t iLightSpot; fListLightSpots->Clear(); for(Int_t iHit=0; iHitGetEntries(); iHit++) { hit = (CbmRichHit*)fListRICHhits->At(iHit); // RICH hit if (strstr(fDebug,"hit")) hit->Print("ugly"); if (hit->GetRefIndex() < 0) continue; // fake hit CbmRichPoint* point = (CbmRichPoint*) fListRICHpts->At(hit->GetRefIndex()); CbmMCTrack* hitTrack = (CbmMCTrack*) fListStack->At(point->GetTrackID()); if(hitTrack->GetPdgCode() != 50000050) continue; // select only Cherenkov photons if (fEbyEgraph) { // fill arrays for hit patterns fXHits->AddAt(hit->X(),fNHits); fYHits->AddAt(hit->Y(),fNHits); } fhPhotonsXY->Fill(hit->X(),hit->Y()); fhNphotPMT ->Fill(hit->GetNPhotons()); fhAmplitude->Fill(hit->GetAmplitude()); fNHits++; // Classify hits by their parents and fill array with light spots Int_t motherTrack = hitTrack->GetMotherID(); Bool_t dejaVue = kFALSE; CbmRichLightSpot *spot; for (iLightSpot=0; iLightSpotAt(iLightSpot); if (spot->GetTrackID() == motherTrack) { dejaVue = kTRUE; spot->AddHitIndex(iHit); break; } } if (!dejaVue) { spot = new((*fListLightSpots)[nLightSpots]) CbmRichLightSpot(); spot->SetTrackID(motherTrack); spot->AddHitIndex(iHit); nLightSpots++; } } if (fNHits > 0) fhNpmt ->Fill(fNHits); // Analyze light spots CbmRichLightSpot *lightSpot; if (nLightSpots > 0) fhNTracks ->Fill(nLightSpots); Int_t nPrimTracks = 0; for (iLightSpot=0; iLightSpotAt(iLightSpot); lightSpot->RingFit(); if (strstr(fDebug,"spot")) lightSpot->Print(); Int_t trackNum = lightSpot->GetTrackID(); CbmMCTrack *cherenkovMother = (CbmMCTrack*) fListStack->At(trackNum); if (cherenkovMother->GetMotherID() <= fNprimary) nPrimTracks++; if (strstr(fDebug,"parent")) { printf("\nCherenkov's mother %d\n",iLightSpot); CbmMCTrack *track; while (1) { printf("==========> Track %d\n",trackNum); track = (CbmMCTrack*)fListStack->At(trackNum); trackNum = track->GetMotherID(); if (trackNum == -1) break; } } fhTrackZ ->Fill(cherenkovMother->GetStartVertex().Z()); Int_t pdg = cherenkovMother->GetPdgCode(); TLorentzVector momentum = cherenkovMother->Get4Momentum(); Double_t pt = momentum.Pt(); Double_t eta = momentum.Eta(); Double_t pabs= momentum.P(); if (TMath::Abs(pdg) == 11) { fhPtElectron ->Fill(pt); fhEtaElectron->Fill(eta); fhPElectron ->Fill(pabs); fhNLSElectron->Fill(pabs,lightSpot->GetNhits()); fhRLSElectron->Fill(pabs,lightSpot->GetRadius()); } else if (TMath::Abs(pdg) == 13) { fhPtMuon ->Fill(pt); fhEtaMuon ->Fill(eta); fhPMuon ->Fill(pabs); fhNLSMuon ->Fill(pabs,lightSpot->GetNhits()); fhRLSMuon ->Fill(pabs,lightSpot->GetRadius()); } else if (TMath::Abs(pdg) == 211) { fhPtPion ->Fill(pt); fhEtaPion ->Fill(eta); fhPPion ->Fill(pabs); fhNLSPion ->Fill(pabs,lightSpot->GetNhits()); fhRLSPion ->Fill(pabs,lightSpot->GetRadius()); } else if (TMath::Abs(pdg) == 321) { fhPtKaon ->Fill(pt); fhEtaKaon ->Fill(eta); fhPKaon ->Fill(pabs); fhNLSKaon ->Fill(pabs,lightSpot->GetNhits()); fhRLSKaon ->Fill(pabs,lightSpot->GetRadius()); } else if (TMath::Abs(pdg) == 2212) { fhPtProton ->Fill(pt); fhEtaProton ->Fill(eta); fhPProton ->Fill(pabs); fhNLSProton ->Fill(pabs,lightSpot->GetNhits()); fhRLSProton ->Fill(pabs,lightSpot->GetRadius()); } } if (nLightSpots > 0) fhNPrimTracks ->Fill(nPrimTracks); // Plot hits per event into a graph and write it to the file if (fEbyEgraph && fNHits > 0) { TGraph *hits = new TGraph(fNHits,fXHits->GetArray(),fYHits->GetArray()); TString name="hitXY_"; name+=fEvent; TString title="Hit pattern in event "; title+=fEvent; hits->SetName (name.Data()); hits->SetTitle(title.Data()); hits->SetMarkerStyle(20); hits->SetMarkerSize(0.2); hits->SetMarkerColor(kBlue); hits->Write(); } } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::AnalyzeRingGuidances() { // Ring guidances analysis if (strstr(fDebug,"ringuid")) printf("\tNumber of RICH ring guidances: %d\n",fListRICHrings->GetEntries()); CbmRichRing* ring=NULL; // Fill arrays with ring centers if (fEbyEgraph) { fNrgs = 0; for (Int_t iRG=0; iRGGetEntries(); iRG++) { ring = (CbmRichRing*)fListRICHrings->At(iRG); // RICH ring guidance // fill arrays for hit patterns fXrgs->AddAt(ring->GetCenterX(),fNrgs); fYrgs->AddAt(ring->GetCenterY(),fNrgs); fNrgs++; } } // Double loop over RICH hits and ring center guidances CbmRichLightSpot *spot; Int_t nLightSpots = fListLightSpots->GetEntries(); for (Int_t iLightSpot=0; iLightSpotAt(iLightSpot); if (spot->GetNhits() < 20) continue; Double_t xSpot = spot->GetCenterX(); Double_t ySpot = spot->GetCenterY(); Double_t rMin = 99.; Int_t idRingMin = 0; Int_t idTrackMin = 0; for (Int_t iRG=0; iRGGetEntries(); iRG++) { ring = (CbmRichRing*)fListRICHrings->At(iRG); Double_t xRing = ring->GetCenterX(); Double_t yRing = ring->GetCenterY(); Float_t rSpotRing = TMath::Sqrt( (xSpot-xRing)*(xSpot-xRing) + (ySpot-yRing)*(ySpot-yRing) ); // search for ring-track combination with shortest distance if (rSpotRing < rMin){ rMin = rSpotRing; idRingMin = spot->GetTrackID(); idTrackMin = ring->GetTrackID(); } if (spot->GetTrackID() == ring->GetTrackID()) { fhRtrueMatch ->Fill(rSpotRing); //pdf(delta_r) for electrons Int_t pdg = ((CbmMCTrack*)fListStack->At(ring->GetTrackID()))->GetPdgCode(); if(pdg == 11 || pdg == -11){ fhPdfElectron->Fill(rSpotRing); } //pdf(delta_r) for pions if(pdg == 211 || pdg == -211){ fhPdfPi->Fill(rSpotRing); } } else fhRfalseMatch->Fill(rSpotRing); } // loop over tracks (ring center guidances) if (rMin != 99.){ if (idRingMin == idTrackMin){ // correct match! fhTrue ->Fill(rMin); Int_t pdg = ((CbmMCTrack*)fListStack->At(idTrackMin))->GetPdgCode(); if (pdg == 11 || pdg == -11) fhTrueElectron ->Fill(rMin); if (pdg == 211 || pdg == -211) fhTruePion ->Fill(rMin); } else { fhFalse ->Fill(rMin); Int_t pdgRing = ((CbmMCTrack*)fListStack->At(idRingMin))->GetPdgCode(); Int_t pdgTrack = ((CbmMCTrack*)fListStack->At(idTrackMin))->GetPdgCode(); if ( TMath::Abs(pdgRing) == 11 && TMath::Abs(pdgTrack) == 211 ) fhFalseElectronPi ->Fill(rMin); if ( TMath::Abs(pdgRing) == 11 ) fhFalseElectron ->Fill(rMin); } } // plot values for closest distance match } // loop over rings (LightSpots) // Plot ring guidances per event into a graph and write it to the file if (fEbyEgraph && fNrgs > 0) { TGraph *rgs = new TGraph(fNrgs,fXrgs->GetArray(),fYrgs->GetArray()); TString name="rgXY_"; name+=fEvent; TString title="Ring guidance pattern in event "; title+=fEvent; rgs->SetName (name.Data()); rgs->SetTitle(title.Data()); rgs->SetMarkerStyle(2); rgs->SetMarkerColor(kRed); rgs->Write(); } fNrgstr = 0; fNrgscp = 0; CbmRichRing* ringtr=NULL; CbmRichRing* ringc=NULL; Double_t x1,y1; if (fEbyEgraph) { for (Int_t iRG2=0; iRG2GetEntries(); iRG2++) { //Points produced by CbmRichProjectionProducer(Stolpovsky) ringc = (CbmRichRing*)fProjectionCollection->At(iRG2); // fill arrays for hit patterns x1=ringc->GetCenterX(); y1=ringc->GetCenterY(); if(x1>-158.&&x1<158.){ if((y1>75.&&y1<250.)||(y1<-75.&&y1>-250.)){ fXcp->AddAt(x1,fNrgscp); fYcp->AddAt(y1,fNrgscp); fNrgscp++; } } } for (Int_t iRG3=0; iRG3GetEntries(); iRG3++) { //Points produced by CbmRichTrackExtrapolation(Hoehne) ringtr = (CbmRichRing*)fTrackCollection->At(iRG3); // fill arrays for hit patterns x1=ringtr->GetCenterX(); y1=ringtr->GetCenterY(); if(x1>-158.&&x1<158.){ if((y1>75.&&y1<250.)||(y1<-75.&&y1>-250.)){ fXtr->AddAt(x1,fNrgstr); fYtr->AddAt(y1,fNrgstr); fNrgstr++; } } } } if (fEbyEgraph && fNrgscp > 0) { TGraph *rgs = new TGraph(fNrgscp,fXcp->GetArray(),fYcp->GetArray()); TString name="cpXY_"; name+=fEvent; TString title="cProjection "; title+=fEvent; rgs->SetName (name.Data()); rgs->SetTitle(title.Data()); rgs->SetMarkerStyle(2); rgs->SetMarkerColor(kRed); rgs->SetMarkerSize(1); rgs->Write(); } if (fEbyEgraph && fNrgstr > 0) { TGraph *rgs = new TGraph(fNrgstr,fXtr->GetArray(),fYtr->GetArray()); TString name="trXY_"; name+=fEvent; TString title="cExtrapolation "; title+=fEvent; rgs->SetName (name.Data()); rgs->SetTitle(title.Data()); rgs->SetMarkerStyle(2); rgs->SetMarkerColor(kGreen); rgs->SetMarkerSize(1); rgs->Write(); } } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::AnalyzeRings() { //Analysis of ring radius resolution. //Convert energy to wavelength Double_t c = 2.998E8; Double_t h = 6.626E-34; Double_t e = 1.6022E-19; Double_t nrefrac = 1.; CbmRichLightSpot* spot; CbmMCPoint* pt; CbmMCTrack* motherTrack; Int_t motherID; Double_t rad, radMC; Double_t lambda_min = 100.; Double_t lambda_max = 700.; for(Int_t iSpot=0; iSpotGetEntries(); iSpot++) { spot = (CbmRichLightSpot*)fListLightSpots->At(iSpot); rad = spot->GetRadius(); motherID = spot->GetTrackID(); motherTrack = (CbmMCTrack*) fListStack->At(motherID); //MC-LightSpot: LightSpot made of MCPoints CbmRichLightSpotMC spotMC; CbmMCTrack* track; Double_t etot,lambda; Int_t nmcp=0; for(Int_t iPoint=0; iPointGetEntries(); iPoint++) { pt = (CbmMCPoint*)fListRICHpts->At(iPoint); TVector3 mom; pt->Momentum(mom); etot = sqrt(mom.Px()*mom.Px()+ mom.Py()*mom.Py() +mom.Pz()*mom.Pz()); lambda = c/nrefrac*h/e/etot; // wavelength in nm (energy in GeV) if(strstr(fDebug,"ring+")) { printf("\tlambda: %f\n", lambda);} track = (CbmMCTrack*) fListStack->At(pt->GetTrackID()); if(track->GetMotherID() == motherID) { nmcp++; if(lambda >= lambda_min && lambda < lambda_max) spotMC.AddMCPointIndex(iPoint); } } if(strstr(fDebug,"ring")) printf("\tNumber of MC point in MC-LightSpot: %d\n", spotMC.GetNPoints()); spotMC.RingFit(); radMC = spotMC.GetRadius(); //Difference between "MCPoint" radius and "fitted". if(strstr(fDebug,"ring")) printf("\tRadius: %f cm. MC radius: %f cm. \n", rad,radMC); if( motherTrack->GetMotherID() == -1 && nmcp>150) //only primary tracks,good rings.. { // fhRadMCXY->Fill(spotMC.GetCenterX(),spotMC.GetCenterY(),radMC); // fhRadHitsXY->Fill(spot->GetCenterX(),spot->GetCenterY(),rad); fhNphotMCRing->Fill(spotMC.GetNPoints()); } } } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::AnalyzeAcceptanceElectrons() { // Acceptance analysis if (!fListStack->GetEntries()) { Fatal("AnalyzeAcceptances","No particles in stack!"); return; } if (((CbmMCTrack*)fListStack->At(0))->GetPdgCode() != 11) return; // if we have the single primary charged particle in event CbmMCTrack* primary = (CbmMCTrack*)fListStack->At(0); Double_t E = primary->Get4Momentum().Energy(); Double_t pz = primary->GetMomentum().Z(); Double_t y = 0.5*TMath::Log( (E+pz)/(E-pz) ); Double_t pt = primary->GetMomentum().Pt(); fhPtYrap->Fill(pt, y); // loop over LightSpots CbmRichLightSpot *spot; Int_t nLightSpots = fListLightSpots->GetEntries(); for (Int_t iLightSpot=0; iLightSpotAt(iLightSpot); if( spot->GetTrackID() != 0 ) continue; if ( spot->GetNhits() < 3 ) continue; // at least 3 hits in ring if ( primary->GetStsPoints() < 3 ) continue; // at least 3 points in STS fhPtYrapSignal->Fill(pt, y); } } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::AnalyzeAcceptanceVectorMesons() { // Acceptance analysis for rho --> e+e- if (!fListStack->GetEntries()) { Fatal("AnalyzeAcceptances","No particles in stack!"); return; } if ( ((CbmMCTrack*)fListStack->At(0))->GetPdgCode() != 113) return; // if we have the single primary charged particle in event CbmMCTrack* primary = (CbmMCTrack*)fListStack->At(0); Double_t E = primary->Get4Momentum().Energy(); Double_t pz = primary->GetMomentum().Z(); Double_t y = 0.5*TMath::Log( (E+pz)/(E-pz) ); Double_t pt = primary->GetMomentum().Pt(); fhPtYrap->Fill(pt, y); Int_t electron1=0; Int_t electron2=0; Float_t rSpotRing; Double_t xSpot, ySpot, xRing, yRing; Int_t iLightSpot,iRG; CbmRichLightSpot *spot; CbmRichRing* ring; Int_t nLightSpots = fListLightSpots->GetEntries(); //electron 1 for ( iLightSpot=0; iLightSpotAt(iLightSpot); if( spot->GetTrackID() != 1 ) continue; if ( (spot->GetNhits() < 17) || (spot->GetNhits() > 40) ) continue; if ( (spot->GetRadius() < 4) || (spot->GetRadius() > 6) ) continue; xSpot = spot->GetCenterX(); ySpot = spot->GetCenterY(); for ( iRG=0; iRGGetEntries(); iRG++) { ring = (CbmRichRing*)fListRICHrings->At(iRG); if(ring->GetTrackID() != 1) continue; xRing = ring->GetCenterX(); yRing = ring->GetCenterY(); rSpotRing = TMath::Sqrt((xSpot-xRing)*(xSpot-xRing) + (ySpot-yRing)*(ySpot-yRing)); if(rSpotRing > 0.77) continue; electron1=1; } } //electron 2 for ( iLightSpot=0; iLightSpotAt(iLightSpot); if( spot->GetTrackID() != 2 ) continue; if ( (spot->GetNhits() < 17) || (spot->GetNhits() > 40) ) continue; if ( (spot->GetRadius() < 4) || (spot->GetRadius() > 6) ) continue; xSpot = spot->GetCenterX(); ySpot = spot->GetCenterY(); for ( iRG=0; iRGGetEntries(); iRG++) { ring = (CbmRichRing*)fListRICHrings->At(iRG); if(ring->GetTrackID() != 2) continue; xRing = ring->GetCenterX(); yRing = ring->GetCenterY(); rSpotRing = TMath::Sqrt((xSpot-xRing)*(xSpot-xRing) + (ySpot-yRing)*(ySpot-yRing)); if(rSpotRing > 0.77) continue; electron2=1; } } if( electron1 && electron2 ) { fhPtYrapSignal->Fill(pt, y); TLorentzVector p1,p2,p; p1 = ((CbmMCTrack*)fListStack->At(0))->Get4Momentum(); p2 = ((CbmMCTrack*)fListStack->At(1))->Get4Momentum(); p = p1 + p2; fhDileptonMass->Fill(p.M()); // inv. mass } } //------------------------------------------------------------------------------ void CbmRichAnalysisHits::Finish() { } //------------------------------------------------------------------------------ ClassImp(CbmRichAnalysisHits)