]>
Commit | Line | Data |
---|---|---|
dc740de4 | 1 | /* $Id$ */ |
2 | ||
79ab56b9 | 3 | #include "AlidNdEtaCorrectionSelector.h" |
539b6cb4 | 4 | |
5 | #include <TStyle.h> | |
6 | #include <TSystem.h> | |
7 | #include <TCanvas.h> | |
8 | #include <TParticle.h> | |
9 | #include <TParticlePDG.h> | |
ad43cff6 | 10 | #include <TH1F.h> |
539b6cb4 | 11 | |
745e836a | 12 | #include <TChain.h> |
13 | #include <TSelector.h> | |
45e97e28 | 14 | #include <TFile.h> |
745e836a | 15 | |
539b6cb4 | 16 | #include <AliLog.h> |
539b6cb4 | 17 | #include <AliTracker.h> |
fcf2fb36 | 18 | #include <AliESDVertex.h> |
19 | #include <AliESD.h> | |
20 | #include <AliESDtrack.h> | |
45e97e28 | 21 | #include <AliRunLoader.h> |
22 | #include <AliStack.h> | |
539b6cb4 | 23 | |
7584d357 | 24 | #include <AliHeader.h> |
25 | #include <AliGenEventHeader.h> | |
4c351225 | 26 | #include <../PYTHIA6/AliGenPythiaEventHeader.h> |
27 | #include <../EVGEN/AliGenCocktailEventHeader.h> | |
7584d357 | 28 | |
37dbb69e | 29 | #include "esdTrackCuts/AliESDtrackCuts.h" |
45e97e28 | 30 | #include "dNdEta/AlidNdEtaCorrection.h" |
31 | #include "AliPWG0Helper.h" | |
4c351225 | 32 | #include "AliPWG0depHelper.h" |
539b6cb4 | 33 | |
74fd10b3 | 34 | #include "dNdEta/dNdEtaAnalysis.h" |
35 | ||
79ab56b9 | 36 | ClassImp(AlidNdEtaCorrectionSelector) |
539b6cb4 | 37 | |
fcf2fb36 | 38 | AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector() : |
16e24ca3 | 39 | AliSelectorRL(), |
539b6cb4 | 40 | fEsdTrackCuts(0), |
406eb63e | 41 | fdNdEtaCorrection(0), |
74fd10b3 | 42 | fdNdEtaAnalysisMC(0), |
43 | fdNdEtaAnalysisESD(0), | |
ad43cff6 | 44 | fPIDParticles(0), |
45 | fPIDTracks(0), | |
61385583 | 46 | fClustersITSPos(0), |
47 | fClustersTPCPos(0), | |
48 | fClustersITSNeg(0), | |
49 | fClustersTPCNeg(0), | |
406eb63e | 50 | fSignMode(0) |
539b6cb4 | 51 | { |
79ab56b9 | 52 | // |
539b6cb4 | 53 | // Constructor. Initialization of pointers |
79ab56b9 | 54 | // |
539b6cb4 | 55 | } |
56 | ||
79ab56b9 | 57 | AlidNdEtaCorrectionSelector::~AlidNdEtaCorrectionSelector() |
539b6cb4 | 58 | { |
79ab56b9 | 59 | // |
60 | // Destructor | |
61 | // | |
539b6cb4 | 62 | |
63 | // histograms are in the output list and deleted when the output | |
64 | // list is deleted by the TSelector dtor | |
65 | } | |
66 | ||
ad43cff6 | 67 | Bool_t AlidNdEtaCorrectionSelector::SignOK(TParticlePDG* particle) |
406eb63e | 68 | { |
69 | // returns if a particle with this sign should be counted | |
70 | // this is determined by the value of fSignMode, which should have the same sign | |
71 | // as the charge | |
72 | // if fSignMode is 0 all particles are counted | |
73 | ||
ad43cff6 | 74 | if (fSignMode == 0) |
75 | return kTRUE; | |
76 | ||
77 | if (!particle) | |
78 | { | |
79 | printf("WARNING: not counting a particle that does not have a pdg particle\n"); | |
80 | return kFALSE; | |
81 | } | |
82 | ||
83 | Double_t charge = particle->Charge(); | |
84 | ||
406eb63e | 85 | if (fSignMode > 0) |
86 | if (charge < 0) | |
87 | return kFALSE; | |
88 | ||
89 | if (fSignMode < 0) | |
90 | if (charge > 0) | |
91 | return kFALSE; | |
92 | ||
93 | return kTRUE; | |
94 | } | |
95 | ||
0bd1f8a0 | 96 | void AlidNdEtaCorrectionSelector::ReadUserObjects(TTree* tree) |
97 | { | |
98 | // read the user objects, called from slavebegin and begin | |
99 | ||
100 | if (!fEsdTrackCuts && fInput) | |
101 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts")); | |
102 | ||
103 | if (!fEsdTrackCuts && tree) | |
104 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts")); | |
105 | ||
106 | if (!fEsdTrackCuts) | |
107 | AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list."); | |
108 | } | |
109 | ||
79ab56b9 | 110 | void AlidNdEtaCorrectionSelector::Begin(TTree * tree) |
539b6cb4 | 111 | { |
112 | // The Begin() function is called at the start of the query. | |
113 | // When running with PROOF Begin() is only called on the client. | |
114 | // The tree argument is deprecated (on PROOF 0 is passed). | |
115 | ||
16e24ca3 | 116 | AliSelectorRL::Begin(tree); |
406eb63e | 117 | |
0bd1f8a0 | 118 | ReadUserObjects(tree); |
119 | ||
406eb63e | 120 | TString option = GetOption(); |
121 | AliInfo(Form("Called with option %s.", option.Data())); | |
122 | ||
123 | if (option.Contains("only-positive")) | |
124 | { | |
125 | AliInfo("Processing only positive particles."); | |
126 | fSignMode = 1; | |
127 | } | |
128 | else if (option.Contains("only-negative")) | |
129 | { | |
130 | AliInfo("Processing only negative particles."); | |
131 | fSignMode = -1; | |
132 | } | |
539b6cb4 | 133 | } |
134 | ||
79ab56b9 | 135 | void AlidNdEtaCorrectionSelector::SlaveBegin(TTree * tree) |
539b6cb4 | 136 | { |
137 | // The SlaveBegin() function is called after the Begin() function. | |
138 | // When running with PROOF SlaveBegin() is called on each slave server. | |
139 | // The tree argument is deprecated (on PROOF 0 is passed). | |
140 | ||
16e24ca3 | 141 | AliSelectorRL::SlaveBegin(tree); |
539b6cb4 | 142 | |
0bd1f8a0 | 143 | ReadUserObjects(tree); |
539b6cb4 | 144 | |
0bd1f8a0 | 145 | fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction"); |
ad43cff6 | 146 | |
147 | fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5); | |
148 | fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5); | |
149 | ||
150 | fClustersITSPos = new TH1F("clusters_its_pos", "clusters_its_pos", 7, -0.5, 6.5); | |
151 | fClustersTPCPos = new TH1F("clusters_tpc_pos", "clusters_tpc_pos", 160, -0.5, 159.5); | |
152 | ||
153 | fClustersITSNeg = new TH1F("clusters_its_neg", "clusters_its_neg", 7, -0.5, 6.5); | |
154 | fClustersTPCNeg = new TH1F("clusters_tpc_neg", "clusters_tpc_neg", 160, -0.5, 159.5); | |
74fd10b3 | 155 | |
156 | fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC"); | |
157 | fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD"); | |
539b6cb4 | 158 | } |
159 | ||
4c351225 | 160 | void AlidNdEtaCorrectionSelector::Init(TTree* tree) |
161 | { | |
162 | // read the user objects | |
163 | ||
164 | AliSelectorRL::Init(tree); | |
165 | ||
166 | // Enable only the needed branches | |
167 | if (tree) | |
168 | { | |
169 | tree->SetBranchStatus("*", 0); | |
170 | tree->SetBranchStatus("fTriggerMask", 1); | |
171 | tree->SetBranchStatus("fSPDVertex*", 1); | |
172 | tree->SetBranchStatus("fTracks.fLabel", 1); | |
173 | tree->SetBranchStatus("fTracks.fITSncls", 1); | |
174 | tree->SetBranchStatus("fTracks.fTPCncls", 1); | |
175 | ||
176 | AliESDtrackCuts::EnableNeededBranches(tree); | |
177 | } | |
178 | } | |
179 | ||
79ab56b9 | 180 | Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry) |
539b6cb4 | 181 | { |
182 | // The Process() function is called for each entry in the tree (or possibly | |
183 | // keyed object in the case of PROOF) to be processed. The entry argument | |
184 | // specifies which entry in the currently loaded tree is to be processed. | |
185 | // It can be passed to either TTree::GetEntry() or TBranch::GetEntry() | |
186 | // to read either all or the required parts of the data. When processing | |
187 | // keyed objects with PROOF, the object is already loaded and is available | |
188 | // via the fObject pointer. | |
189 | // | |
190 | // This function should contain the "body" of the analysis. It can contain | |
191 | // simple or elaborate selection criteria, run algorithms on the data | |
192 | // of the event and typically fill histograms. | |
193 | ||
194 | // WARNING when a selector is used with a TChain, you must use | |
195 | // the pointer to the current TTree to call GetEntry(entry). | |
196 | // The entry is always the local entry number in the current tree. | |
16e24ca3 | 197 | // Assuming that fTree is the pointer to the TChain being processed, |
198 | // use fTree->GetTree()->GetEntry(entry). | |
539b6cb4 | 199 | |
7584d357 | 200 | AliDebug(AliLog::kDebug+1,"Processing event ...\n"); |
201 | ||
202 | ||
16e24ca3 | 203 | if (AliSelectorRL::Process(entry) == kFALSE) |
539b6cb4 | 204 | return kFALSE; |
205 | ||
b8e8168f | 206 | // check prerequesites |
207 | if (!fESD) | |
208 | { | |
209 | AliDebug(AliLog::kError, "ESD branch not available"); | |
210 | return kFALSE; | |
211 | } | |
212 | ||
dc740de4 | 213 | AliHeader* header = GetHeader(); |
214 | if (!header) | |
b8e8168f | 215 | { |
dc740de4 | 216 | AliDebug(AliLog::kError, "Header not available"); |
b8e8168f | 217 | return kFALSE; |
218 | } | |
219 | ||
d09fb536 | 220 | AliStack* stack = GetStack(); |
45e97e28 | 221 | if (!stack) |
222 | { | |
223 | AliDebug(AliLog::kError, "Stack not available"); | |
224 | return kFALSE; | |
225 | } | |
226 | ||
b8e8168f | 227 | if (!fEsdTrackCuts) |
228 | { | |
229 | AliDebug(AliLog::kError, "fESDTrackCuts not available"); | |
539b6cb4 | 230 | return kFALSE; |
b8e8168f | 231 | } |
539b6cb4 | 232 | |
45e97e28 | 233 | Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD); |
745e836a | 234 | |
45e97e28 | 235 | Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD); |
539b6cb4 | 236 | |
539b6cb4 | 237 | // get the MC vertex |
dc740de4 | 238 | AliGenEventHeader* genHeader = header->GenEventHeader(); |
539b6cb4 | 239 | |
74fd10b3 | 240 | // primary vertex (from MC) |
539b6cb4 | 241 | TArrayF vtxMC(3); |
242 | genHeader->PrimaryVertex(vtxMC); | |
539b6cb4 | 243 | |
74fd10b3 | 244 | // get process type |
245 | Int_t processType = AliPWG0depHelper::GetPythiaEventProcessType(header); | |
246 | AliDebug(AliLog::kDebug+1, Form("Found pythia procces type %d", processType)); | |
247 | ||
248 | if (processType<0) | |
249 | AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType)); | |
250 | ||
539b6cb4 | 251 | // loop over mc particles |
45e97e28 | 252 | Int_t nPrim = stack->GetNprimary(); |
539b6cb4 | 253 | |
45e97e28 | 254 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) |
539b6cb4 | 255 | { |
45e97e28 | 256 | AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc)); |
257 | ||
258 | TParticle* particle = stack->Particle(iMc); | |
539b6cb4 | 259 | |
260 | if (!particle) | |
45e97e28 | 261 | { |
262 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc)); | |
539b6cb4 | 263 | continue; |
45e97e28 | 264 | } |
539b6cb4 | 265 | |
45e97e28 | 266 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) |
539b6cb4 | 267 | continue; |
268 | ||
ad43cff6 | 269 | if (SignOK(particle->GetPDG()) == kFALSE) |
10ebe68d | 270 | continue; |
ad43cff6 | 271 | |
745e836a | 272 | Float_t eta = particle->Eta(); |
45e97e28 | 273 | Float_t pt = particle->Pt(); |
274 | ||
74fd10b3 | 275 | fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, processType); |
7584d357 | 276 | |
74fd10b3 | 277 | if (eventTriggered) |
278 | { | |
279 | if (vertexReconstructed) | |
280 | { | |
281 | fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt); | |
282 | } | |
ad43cff6 | 283 | } |
74fd10b3 | 284 | } // end of mc particle |
e8cb44f1 | 285 | |
539b6cb4 | 286 | // ######################################################## |
287 | // loop over esd tracks | |
288 | Int_t nTracks = fESD->GetNumberOfTracks(); | |
45e97e28 | 289 | |
406eb63e | 290 | // count the number of "good" tracks as parameter for vertex reconstruction efficiency |
45e97e28 | 291 | Int_t nGoodTracks = 0; |
539b6cb4 | 292 | for (Int_t t=0; t<nTracks; t++) |
293 | { | |
45e97e28 | 294 | AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t)); |
295 | ||
539b6cb4 | 296 | AliESDtrack* esdTrack = fESD->GetTrack(t); |
297 | ||
298 | // cut the esd track? | |
299 | if (!fEsdTrackCuts->AcceptTrack(esdTrack)) | |
300 | continue; | |
301 | ||
45e97e28 | 302 | nGoodTracks++; |
303 | ||
304 | // using the properties of the mc particle | |
539b6cb4 | 305 | Int_t label = TMath::Abs(esdTrack->GetLabel()); |
79ab56b9 | 306 | if (label == 0) |
539b6cb4 | 307 | { |
b8e8168f | 308 | AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t)); |
539b6cb4 | 309 | continue; |
310 | } | |
539b6cb4 | 311 | |
45e97e28 | 312 | TParticle* particle = stack->Particle(label); |
313 | if (!particle) | |
314 | { | |
315 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label)); | |
316 | continue; | |
317 | } | |
539b6cb4 | 318 | |
ad43cff6 | 319 | if (SignOK(particle->GetPDG()) == kFALSE) |
406eb63e | 320 | continue; |
321 | ||
74fd10b3 | 322 | if (eventTriggered && vertexReconstructed) |
ad43cff6 | 323 | { |
74fd10b3 | 324 | fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt()); |
325 | fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt()); | |
ad43cff6 | 326 | if (particle->Pt() > 0.1 && particle->Pt() < 0.2) |
74fd10b3 | 327 | { |
328 | fPIDTracks->Fill(particle->GetPdgCode()); | |
ad43cff6 | 329 | if (particle->GetPDG()->Charge() > 0) |
330 | { | |
331 | fClustersITSPos->Fill(esdTrack->GetITSclusters(0)); | |
332 | fClustersTPCPos->Fill(esdTrack->GetTPCclusters(0)); | |
333 | } | |
334 | else | |
335 | { | |
336 | fClustersITSNeg->Fill(esdTrack->GetITSclusters(0)); | |
337 | fClustersTPCNeg->Fill(esdTrack->GetTPCclusters(0)); | |
338 | } | |
339 | } | |
340 | } | |
539b6cb4 | 341 | } // end of track loop |
342 | ||
74fd10b3 | 343 | if (eventTriggered && vertexReconstructed) |
344 | fdNdEtaAnalysisMC->FillEvent(vtxMC[2], nGoodTracks); | |
345 | ||
7584d357 | 346 | // stuff regarding the vertex reco correction and trigger bias correction |
74fd10b3 | 347 | fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, processType); |
7584d357 | 348 | if (eventTriggered) { |
847489f7 | 349 | if (vertexReconstructed) |
74fd10b3 | 350 | { |
351 | fdNdEtaAnalysisESD->FillEvent(vtxMC[2], nGoodTracks); | |
352 | } | |
847489f7 | 353 | } |
45e97e28 | 354 | |
539b6cb4 | 355 | return kTRUE; |
356 | } | |
357 | ||
79ab56b9 | 358 | void AlidNdEtaCorrectionSelector::SlaveTerminate() |
539b6cb4 | 359 | { |
360 | // The SlaveTerminate() function is called after all entries or objects | |
361 | // have been processed. When running with PROOF SlaveTerminate() is called | |
362 | // on each slave server. | |
363 | ||
16e24ca3 | 364 | AliSelectorRL::SlaveTerminate(); |
539b6cb4 | 365 | |
366 | // Add the histograms to the output on each slave server | |
367 | if (!fOutput) | |
368 | { | |
b8e8168f | 369 | AliDebug(AliLog::kError, "ERROR: Output list not initialized"); |
539b6cb4 | 370 | return; |
371 | } | |
79ab56b9 | 372 | |
745e836a | 373 | fOutput->Add(fdNdEtaCorrection); |
74fd10b3 | 374 | fOutput->Add(fdNdEtaAnalysisMC); |
375 | fOutput->Add(fdNdEtaAnalysisESD); | |
539b6cb4 | 376 | } |
377 | ||
79ab56b9 | 378 | void AlidNdEtaCorrectionSelector::Terminate() |
539b6cb4 | 379 | { |
380 | // The Terminate() function is the last function to be called during | |
381 | // a query. It always runs on the client, it can be used to present | |
382 | // the results graphically or save the results to file. | |
383 | ||
16e24ca3 | 384 | AliSelectorRL::Terminate(); |
539b6cb4 | 385 | |
45e97e28 | 386 | fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction")); |
74fd10b3 | 387 | fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC")); |
388 | fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD")); | |
389 | if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD) | |
ad43cff6 | 390 | { |
391 | AliDebug(AliLog::kError, "Could not read object from output list"); | |
392 | return; | |
393 | } | |
79ab56b9 | 394 | |
45e97e28 | 395 | fdNdEtaCorrection->Finish(); |
539b6cb4 | 396 | |
406eb63e | 397 | TFile* fout = new TFile(Form("correction_map%s.root", GetOption()), "RECREATE"); |
539b6cb4 | 398 | |
0bd1f8a0 | 399 | if (fEsdTrackCuts) |
400 | fEsdTrackCuts->SaveHistograms("esd_track_cuts"); | |
45e97e28 | 401 | fdNdEtaCorrection->SaveHistograms(); |
74fd10b3 | 402 | fdNdEtaAnalysisMC->SaveHistograms(); |
403 | fdNdEtaAnalysisESD->SaveHistograms(); | |
539b6cb4 | 404 | |
405 | fout->Write(); | |
406 | fout->Close(); | |
407 | ||
45e97e28 | 408 | fdNdEtaCorrection->DrawHistograms(); |
ad43cff6 | 409 | |
0bd1f8a0 | 410 | if (fPIDParticles && fPIDTracks) |
411 | { | |
412 | new TCanvas("pidcanvas", "pidcanvas", 500, 500); | |
ad43cff6 | 413 | |
0bd1f8a0 | 414 | fPIDParticles->Draw(); |
415 | fPIDTracks->SetLineColor(2); | |
416 | fPIDTracks->Draw("SAME"); | |
ad43cff6 | 417 | |
0bd1f8a0 | 418 | TDatabasePDG* pdgDB = new TDatabasePDG; |
ad43cff6 | 419 | |
0bd1f8a0 | 420 | for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i) |
421 | if (fPIDParticles->GetBinContent(i) > 0) | |
422 | 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)); | |
ad43cff6 | 423 | |
0bd1f8a0 | 424 | delete pdgDB; |
425 | pdgDB = 0; | |
426 | } | |
ad43cff6 | 427 | |
0bd1f8a0 | 428 | if (fClustersITSPos && fClustersITSNeg && fClustersTPCPos && fClustersTPCNeg) |
429 | { | |
430 | TCanvas* canvas = new TCanvas("clusters", "clusters", 1000, 500); | |
431 | canvas->Divide(2, 1); | |
432 | ||
433 | canvas->cd(1); | |
434 | fClustersITSPos->Draw(); | |
435 | fClustersITSNeg->SetLineColor(kRed); | |
436 | fClustersITSNeg->Draw("SAME"); | |
437 | ||
438 | canvas->cd(2); | |
439 | fClustersTPCPos->Draw(); | |
440 | fClustersTPCNeg->SetLineColor(kRed); | |
441 | fClustersTPCNeg->Draw("SAME"); | |
442 | } | |
539b6cb4 | 443 | } |