3 #include "AlidNdEtaCorrectionSelector.h"
9 #include <TParticlePDG.h>
13 #include <TSelector.h>
17 #include <AliGenEventHeader.h>
18 #include <AliTracker.h>
19 #include <AliHeader.h>
20 #include <AliESDVertex.h>
22 #include <AliESDtrack.h>
23 #include <AliRunLoader.h>
26 #include "esdTrackCuts/AliESDtrackCuts.h"
27 #include "dNdEta/AlidNdEtaCorrection.h"
28 #include "AliPWG0Helper.h"
30 ClassImp(AlidNdEtaCorrectionSelector)
32 AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector() :
45 // Constructor. Initialization of pointers
49 AlidNdEtaCorrectionSelector::~AlidNdEtaCorrectionSelector()
55 // histograms are in the output list and deleted when the output
56 // list is deleted by the TSelector dtor
59 Bool_t AlidNdEtaCorrectionSelector::SignOK(TParticlePDG* particle)
61 // returns if a particle with this sign should be counted
62 // this is determined by the value of fSignMode, which should have the same sign
64 // if fSignMode is 0 all particles are counted
71 printf("WARNING: not counting a particle that does not have a pdg particle\n");
75 Double_t charge = particle->Charge();
88 void AlidNdEtaCorrectionSelector::Begin(TTree * tree)
90 // The Begin() function is called at the start of the query.
91 // When running with PROOF Begin() is only called on the client.
92 // The tree argument is deprecated (on PROOF 0 is passed).
94 AliSelectorRL::Begin(tree);
96 TString option = GetOption();
97 AliInfo(Form("Called with option %s.", option.Data()));
99 if (option.Contains("only-positive"))
101 AliInfo("Processing only positive particles.");
104 else if (option.Contains("only-negative"))
106 AliInfo("Processing only negative particles.");
111 void AlidNdEtaCorrectionSelector::SlaveBegin(TTree * tree)
113 // The SlaveBegin() function is called after the Begin() function.
114 // When running with PROOF SlaveBegin() is called on each slave server.
115 // The tree argument is deprecated (on PROOF 0 is passed).
117 AliSelectorRL::SlaveBegin(tree);
119 fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
122 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fTree->GetUserInfo()->FindObject("AliESDtrackCuts"));
125 AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from user info");
127 fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
128 fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
130 fClustersITSPos = new TH1F("clusters_its_pos", "clusters_its_pos", 7, -0.5, 6.5);
131 fClustersTPCPos = new TH1F("clusters_tpc_pos", "clusters_tpc_pos", 160, -0.5, 159.5);
133 fClustersITSNeg = new TH1F("clusters_its_neg", "clusters_its_neg", 7, -0.5, 6.5);
134 fClustersTPCNeg = new TH1F("clusters_tpc_neg", "clusters_tpc_neg", 160, -0.5, 159.5);
137 Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
139 // The Process() function is called for each entry in the tree (or possibly
140 // keyed object in the case of PROOF) to be processed. The entry argument
141 // specifies which entry in the currently loaded tree is to be processed.
142 // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
143 // to read either all or the required parts of the data. When processing
144 // keyed objects with PROOF, the object is already loaded and is available
145 // via the fObject pointer.
147 // This function should contain the "body" of the analysis. It can contain
148 // simple or elaborate selection criteria, run algorithms on the data
149 // of the event and typically fill histograms.
151 // WARNING when a selector is used with a TChain, you must use
152 // the pointer to the current TTree to call GetEntry(entry).
153 // The entry is always the local entry number in the current tree.
154 // Assuming that fTree is the pointer to the TChain being processed,
155 // use fTree->GetTree()->GetEntry(entry).
157 if (AliSelectorRL::Process(entry) == kFALSE)
160 // check prerequesites
163 AliDebug(AliLog::kError, "ESD branch not available");
167 AliHeader* header = GetHeader();
170 AliDebug(AliLog::kError, "Header not available");
174 AliStack* stack = GetStack();
177 AliDebug(AliLog::kError, "Stack not available");
183 AliDebug(AliLog::kError, "fESDTrackCuts not available");
187 Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
189 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
191 fdNdEtaCorrection->IncreaseEventCount();
193 fdNdEtaCorrection->IncreaseTriggeredEventCount();
196 AliGenEventHeader* genHeader = header->GenEventHeader();
199 genHeader->PrimaryVertex(vtxMC);
201 // loop over mc particles
202 Int_t nPrim = stack->GetNprimary();
204 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
206 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
208 TParticle* particle = stack->Particle(iMc);
212 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
216 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
219 if (SignOK(particle->GetPDG()) == kFALSE)
222 Float_t eta = particle->Eta();
223 Float_t pt = particle->Pt();
225 if (vertexReconstructed)
227 fdNdEtaCorrection->FillParticle(vtxMC[2], eta, pt);
228 if (pt > 0.1 && pt < 0.2)
229 fPIDParticles->Fill(particle->GetPdgCode());
232 fdNdEtaCorrection->FillParticleAllEvents(eta, pt);
234 fdNdEtaCorrection->FillParticleWhenEventTriggered(eta, pt);
235 }// end of mc particle
237 // ########################################################
238 // loop over esd tracks
239 Int_t nTracks = fESD->GetNumberOfTracks();
241 // count the number of "good" tracks as parameter for vertex reconstruction efficiency
242 Int_t nGoodTracks = 0;
243 for (Int_t t=0; t<nTracks; t++)
245 AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t));
247 AliESDtrack* esdTrack = fESD->GetTrack(t);
249 // cut the esd track?
250 if (!fEsdTrackCuts->AcceptTrack(esdTrack))
255 // using the properties of the mc particle
256 Int_t label = TMath::Abs(esdTrack->GetLabel());
259 AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
263 TParticle* particle = stack->Particle(label);
266 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
270 if (SignOK(particle->GetPDG()) == kFALSE)
273 if (vertexReconstructed)
275 fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
276 if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
278 fPIDTracks->Fill(particle->GetPdgCode());
279 if (particle->GetPDG()->Charge() > 0)
281 fClustersITSPos->Fill(esdTrack->GetITSclusters(0));
282 fClustersTPCPos->Fill(esdTrack->GetTPCclusters(0));
286 fClustersITSNeg->Fill(esdTrack->GetITSclusters(0));
287 fClustersTPCNeg->Fill(esdTrack->GetTPCclusters(0));
291 } // end of track loop
293 fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
296 fdNdEtaCorrection->FillEventWithTrigger(vtxMC[2], nGoodTracks);
297 if (vertexReconstructed)
298 fdNdEtaCorrection->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
304 void AlidNdEtaCorrectionSelector::SlaveTerminate()
306 // The SlaveTerminate() function is called after all entries or objects
307 // have been processed. When running with PROOF SlaveTerminate() is called
308 // on each slave server.
310 AliSelectorRL::SlaveTerminate();
312 // Add the histograms to the output on each slave server
315 AliDebug(AliLog::kError, "ERROR: Output list not initialized");
319 fOutput->Add(fdNdEtaCorrection);
322 void AlidNdEtaCorrectionSelector::Terminate()
324 // The Terminate() function is the last function to be called during
325 // a query. It always runs on the client, it can be used to present
326 // the results graphically or save the results to file.
328 AliSelectorRL::Terminate();
330 fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
331 if (!fdNdEtaCorrection)
333 AliDebug(AliLog::kError, "Could not read object from output list");
337 fdNdEtaCorrection->Finish();
339 TFile* fout = new TFile(Form("correction_map%s.root", GetOption()), "RECREATE");
341 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
342 fdNdEtaCorrection->SaveHistograms();
347 fdNdEtaCorrection->DrawHistograms();
349 new TCanvas("pidcanvas", "pidcanvas", 500, 500);
351 fPIDParticles->Draw();
352 fPIDTracks->SetLineColor(2);
353 fPIDTracks->Draw("SAME");
355 TDatabasePDG* pdgDB = new TDatabasePDG;
357 for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
358 if (fPIDParticles->GetBinContent(i) > 0)
359 printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
364 TCanvas* canvas = new TCanvas("clusters", "clusters", 1000, 500);
365 canvas->Divide(2, 1);
368 fClustersITSPos->Draw();
369 fClustersITSNeg->SetLineColor(kRed);
370 fClustersITSNeg->Draw("SAME");
373 fClustersTPCPos->Draw();
374 fClustersTPCNeg->SetLineColor(kRed);
375 fClustersTPCNeg->Draw("SAME");