fixing warnings
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / AliMultiplicityTask.cxx
CommitLineData
a9017e49 1/* $Id$ */
2
3#include "AliMultiplicityTask.h"
4
5#include <TStyle.h>
6#include <TSystem.h>
7#include <TCanvas.h>
8#include <TVector3.h>
9#include <TChain.h>
10#include <TFile.h>
11#include <TH2F.h>
12#include <TH3F.h>
13#include <TParticle.h>
14#include <TRandom.h>
15#include <TNtuple.h>
16#include <TObjString.h>
17#include <TF1.h>
18
19#include <AliLog.h>
20#include <AliESDVertex.h>
21#include <AliESDEvent.h>
22#include <AliStack.h>
23#include <AliHeader.h>
24#include <AliGenEventHeader.h>
25#include <AliMultiplicity.h>
26#include <AliAnalysisManager.h>
27#include <AliMCEventHandler.h>
28#include <AliMCEvent.h>
29#include <AliESDInputHandler.h>
30
745d6088 31#include "AliESDtrackCuts.h"
a9017e49 32#include "AliPWG0Helper.h"
33#include "dNdEta/AliMultiplicityCorrection.h"
34#include "AliCorrection.h"
35#include "AliCorrectionMatrix3D.h"
36
a9017e49 37ClassImp(AliMultiplicityTask)
38
39AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
40 AliAnalysisTask("AliMultiplicityTask", ""),
41 fESD(0),
42 fOption(opt),
745d6088 43 fAnalysisMode(AliPWG0Helper::kSPD),
2fa65f52 44 fReadMC(kFALSE),
a9017e49 45 fMultiplicity(0),
46 fEsdTrackCuts(0),
47 fSystSkipParticles(kFALSE),
48 fSelectProcessType(0),
49 fParticleSpecies(0),
50 fPtSpectrum(0),
51 fOutput(0)
52{
53 //
54 // Constructor. Initialization of pointers
55 //
56
57 for (Int_t i = 0; i<4; i++)
58 fParticleCorrection[i] = 0;
59
60 // Define input and output slots here
61 DefineInput(0, TChain::Class());
62 DefineOutput(0, TList::Class());
63}
64
65AliMultiplicityTask::~AliMultiplicityTask()
66{
67 //
68 // Destructor
69 //
70
71 // histograms are in the output list and deleted when the output
72 // list is deleted by the TSelector dtor
73
74 if (fOutput) {
75 delete fOutput;
76 fOutput = 0;
77 }
78}
79
80//________________________________________________________________________
81void AliMultiplicityTask::ConnectInputData(Option_t *)
82{
83 // Connect ESD
84 // Called once
85
2fa65f52 86 Printf("AliMultiplicityTask::ConnectInputData called");
87
a9017e49 88 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
89 if (!tree) {
90 Printf("ERROR: Could not read tree from input slot 0");
91 } else {
92 // Disable all branches and enable only the needed ones
2fa65f52 93 //tree->SetBranchStatus("*", 0);
94
a9017e49 95 tree->SetBranchStatus("fTriggerMask", 1);
96 tree->SetBranchStatus("fSPDVertex*", 1);
97
745d6088 98 if (fAnalysisMode == AliPWG0Helper::kSPD)
a9017e49 99 tree->SetBranchStatus("fSPDMult*", 1);
a9017e49 100
745d6088 101 if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
a9017e49 102 AliESDtrackCuts::EnableNeededBranches(tree);
103 tree->SetBranchStatus("fTracks.fLabel", 1);
2fa65f52 104 }
a9017e49 105
a9017e49 106 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
107
108 if (!esdH) {
109 Printf("ERROR: Could not get ESDInputHandler");
110 } else
111 fESD = esdH->GetEvent();
112 }
113}
114
115void AliMultiplicityTask::CreateOutputObjects()
116{
117 // create result objects and add to output list
118
119 fOutput = new TList;
120 fOutput->SetOwner();
121
122 fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
123 fOutput->Add(fMultiplicity);
124
125 if (fOption.Contains("skip-particles"))
126 {
127 fSystSkipParticles = kTRUE;
128 AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
129 }
130
131 if (fOption.Contains("particle-efficiency"))
132 for (Int_t i = 0; i<4; i++)
133 {
134 fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
135 fOutput->Add(fParticleCorrection[i]);
136 }
137
138 if (fOption.Contains("only-process-type-nd"))
139 fSelectProcessType = 1;
140
141 if (fOption.Contains("only-process-type-sd"))
142 fSelectProcessType = 2;
143
144 if (fOption.Contains("only-process-type-dd"))
145 fSelectProcessType = 3;
146
147 if (fSelectProcessType != 0)
148 AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
149
150 if (fOption.Contains("pt-spectrum-hist"))
151 {
152 TFile* file = TFile::Open("ptspectrum_fit.root");
153 if (file)
154 {
155 TString subStr(fOption(fOption.Index("pt-spectrum")+17, 3));
156 TString histName(Form("ptspectrum_%s", subStr.Data()));
157 AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
158 fPtSpectrum = (TH1*) file->Get(histName);
159 if (!fPtSpectrum)
160 AliError("Histogram not found");
161 }
162 else
163 AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
164
165 if (fPtSpectrum)
166 AliInfo("WARNING: Systematic study enabled. Pt spectrum will be modified");
167 }
168
169 if (fOption.Contains("pt-spectrum-func"))
170 {
171 if (fPtSpectrum)
172 {
173 Printf("Using function from input list for systematic p_t study");
174 }
175 else
176 {
177 fPtSpectrum = new TH1F("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
178 fPtSpectrum->SetBinContent(1, 1);
179 }
180
181 if (fPtSpectrum)
182 AliInfo("WARNING: Systematic study enabled. Pt spectrum will be modified");
183 }
184
185 if (fOption.Contains("particle-species")) {
186 fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "vtx:Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec:nolabel:doublePrim:doubleCount");
187 fOutput->Add(fParticleSpecies);
188 }
189
190 // TODO set seed for random generator
191}
192
193void AliMultiplicityTask::Exec(Option_t*)
194{
195 // process the event
196
197 // Check prerequisites
198 if (!fESD)
199 {
745d6088 200 AliDebug(AliLog::kError, "ESD not available");
a9017e49 201 return;
202 }
203
745d6088 204 //MB1 definition
205 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1);
2fa65f52 206 // only FASTOR
745d6088 207 //Bool_t eventTriggered = fESD->GetTriggerMask() & 32;
a9017e49 208
209 Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD->GetVertex());
210
211 // get the ESD vertex
212 const AliESDVertex* vtxESD = fESD->GetVertex();
213 Double_t vtx[3];
214 vtxESD->GetXYZ(vtx);
215
2fa65f52 216 // post the data already here
217 PostData(0, fOutput);
a9017e49 218
219 //const Float_t kPtCut = 0.3;
220
2fa65f52 221 // create list of (label, eta) tuples
222 Int_t inputCount = 0;
223 Int_t* labelArr = 0;
224 Float_t* etaArr = 0;
745d6088 225 if (fAnalysisMode == AliPWG0Helper::kSPD)
a9017e49 226 {
2fa65f52 227 // get tracklets
228 const AliMultiplicity* mult = fESD->GetMultiplicity();
229 if (!mult)
a9017e49 230 {
2fa65f52 231 AliDebug(AliLog::kError, "AliMultiplicity not available");
232 return;
a9017e49 233 }
234
2fa65f52 235 labelArr = new Int_t[mult->GetNumberOfTracklets()];
236 etaArr = new Float_t[mult->GetNumberOfTracklets()];
237
238 // get multiplicity from ITS tracklets
239 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
a9017e49 240 {
2fa65f52 241 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
a9017e49 242
2fa65f52 243 // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
244 if (mult->GetDeltaPhi(i) < -1000)
245 continue;
a9017e49 246
2fa65f52 247 etaArr[inputCount] = mult->GetEta(i);
3886f106 248 // TODO add second label array
249 labelArr[inputCount] = mult->GetLabel(i, 0);
2fa65f52 250 ++inputCount;
251 }
252 }
745d6088 253 else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
2fa65f52 254 {
255 if (!fEsdTrackCuts)
256 {
257 AliDebug(AliLog::kError, "fESDTrackCuts not available");
258 return;
259 }
a9017e49 260
2fa65f52 261 // get multiplicity from ESD tracks
262 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
263 Int_t nGoodTracks = list->GetEntries();
a9017e49 264
2fa65f52 265 labelArr = new Int_t[nGoodTracks];
266 etaArr = new Float_t[nGoodTracks];
a9017e49 267
2fa65f52 268 // loop over esd tracks
269 for (Int_t i=0; i<nGoodTracks; i++)
a9017e49 270 {
2fa65f52 271 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
272 if (!esdTrack)
a9017e49 273 {
2fa65f52 274 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
275 continue;
a9017e49 276 }
2fa65f52 277
278 etaArr[inputCount] = esdTrack->Eta();
279 labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
280 ++inputCount;
a9017e49 281 }
282
2fa65f52 283 delete list;
284 }
a9017e49 285
2fa65f52 286 // eta range for nMCTracksSpecies, nESDTracksSpecies
287 Float_t etaRange = -1;
288 switch (fAnalysisMode) {
745d6088 289 case AliPWG0Helper::kTPC:
290 case AliPWG0Helper::kTPCITS:
291 etaRange = 0.9; break;
292 case AliPWG0Helper::kSPD: etaRange = 2.0; break;
2fa65f52 293 }
a9017e49 294
2fa65f52 295 if (!fReadMC) // Processing of ESD information
296 {
297 if (eventTriggered && eventVertex)
298 {
299 Int_t nESDTracks05 = 0;
300 Int_t nESDTracks09 = 0;
301 Int_t nESDTracks15 = 0;
302 Int_t nESDTracks20 = 0;
a9017e49 303
2fa65f52 304 for (Int_t i=0; i<inputCount; ++i)
305 {
306 Float_t eta = etaArr[i];
a9017e49 307
2fa65f52 308 if (TMath::Abs(eta) < 0.5)
309 nESDTracks05++;
a9017e49 310
2fa65f52 311 if (TMath::Abs(eta) < 0.9)
312 nESDTracks09++;
a9017e49 313
2fa65f52 314 if (TMath::Abs(eta) < 1.5)
315 nESDTracks15++;
316
317 if (TMath::Abs(eta) < 2.0)
318 nESDTracks20++;
a9017e49 319 }
2fa65f52 320
321 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
322 }
323 }
324 else if (fReadMC) // Processing of MC information
325 {
326 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
327 if (!eventHandler) {
328 Printf("ERROR: Could not retrieve MC event handler");
329 return;
a9017e49 330 }
a9017e49 331
2fa65f52 332 AliMCEvent* mcEvent = eventHandler->MCEvent();
333 if (!mcEvent) {
334 Printf("ERROR: Could not retrieve MC event");
335 return;
336 }
a9017e49 337
2fa65f52 338 AliStack* stack = mcEvent->Stack();
339 if (!stack)
340 {
341 AliDebug(AliLog::kError, "Stack not available");
342 return;
343 }
344
345 AliHeader* header = mcEvent->Header();
346 if (!header)
347 {
348 AliDebug(AliLog::kError, "Header not available");
349 return;
350 }
a9017e49 351
2fa65f52 352 Bool_t processEvent = kTRUE;
353 if (fSelectProcessType > 0)
354 {
355 // getting process information; NB: this only works for Pythia
356 Int_t processtype = AliPWG0Helper::GetPythiaEventProcessType(header);
a9017e49 357
2fa65f52 358 processEvent = kFALSE;
a9017e49 359
2fa65f52 360 // non diffractive
361 if (fSelectProcessType == 1 && processtype !=92 && processtype !=93 && processtype != 94)
362 processEvent = kTRUE;
a9017e49 363
2fa65f52 364 // single diffractive
365 if (fSelectProcessType == 2 && (processtype == 92 || processtype == 93))
366 processEvent = kTRUE;
a9017e49 367
2fa65f52 368 // double diffractive
369 if (fSelectProcessType == 3 && processtype == 94)
370 processEvent = kTRUE;
a9017e49 371
2fa65f52 372 if (!processEvent)
373 AliDebug(AliLog::kDebug, Form("Skipping this event, because it is not of the requested process type (%d)", processtype));
374 }
a9017e49 375
376 // systematic study: 10% lower efficiency
377 if (fSystSkipParticles && (gRandom->Uniform() < 0.1))
2fa65f52 378 processEvent = kFALSE;
a9017e49 379
2fa65f52 380 if (processEvent)
381 {
382 // get the MC vertex
383 AliGenEventHeader* genHeader = header->GenEventHeader();
384 TArrayF vtxMC(3);
385 genHeader->PrimaryVertex(vtxMC);
386
387 // get number of tracks from MC
388 // loop over mc particles
389 Int_t nPrim = stack->GetNprimary();
390 Int_t nMCPart = stack->GetNtrack();
391
392 // tracks in different eta ranges
393 Int_t nMCTracks05 = 0;
394 Int_t nMCTracks09 = 0;
395 Int_t nMCTracks15 = 0;
396 Int_t nMCTracks20 = 0;
397 Int_t nMCTracksAll = 0;
398
399 // tracks per particle species, in |eta| < 2 (systematic study)
400 Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
401 for (Int_t i = 0; i<4; ++i)
402 nMCTracksSpecies[i] = 0;
403
404 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
405 {
406 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
a9017e49 407
2fa65f52 408 TParticle* particle = stack->Particle(iMc);
a9017e49 409
2fa65f52 410 if (!particle)
411 {
412 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
413 continue;
414 }
a9017e49 415
2fa65f52 416 Bool_t debug = kFALSE;
417 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
418 {
419 //printf("%d) DROPPED ", iMC);
420 //particle->Print();
421 continue;
422 }
a9017e49 423
2fa65f52 424 //printf("%d) OK ", iMC);
425 //particle->Print();
a9017e49 426
2fa65f52 427 //if (particle->Pt() < kPtCut)
428 // continue;
a9017e49 429
2fa65f52 430 Int_t particleWeight = 1;
a9017e49 431
2fa65f52 432 Float_t pt = particle->Pt();
a9017e49 433
2fa65f52 434 // in case of systematic study, weight according to the change of the pt spectrum
435 // it cannot be just multiplied because we cannot count "half of a particle"
436 // instead a random generator decides if the particle is counted twice (if value > 1)
437 // or not (if value < 0)
438 if (fPtSpectrum)
a9017e49 439 {
2fa65f52 440 Int_t bin = fPtSpectrum->FindBin(pt);
441 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
a9017e49 442 {
2fa65f52 443 Float_t factor = fPtSpectrum->GetBinContent(bin);
444 if (factor > 0)
445 {
446 Float_t random = gRandom->Uniform();
447 if (factor > 1 && random < factor - 1)
448 {
449 particleWeight = 2;
450 }
451 else if (factor < 1 && random < 1 - factor)
452 particleWeight = 0;
453 }
a9017e49 454 }
a9017e49 455 }
a9017e49 456
2fa65f52 457 //Printf("MC weight is: %d", particleWeight);
a9017e49 458
2fa65f52 459 if (TMath::Abs(particle->Eta()) < 0.5)
460 nMCTracks05 += particleWeight;
a9017e49 461
2fa65f52 462 if (TMath::Abs(particle->Eta()) < 0.9)
463 nMCTracks09 += particleWeight;
a9017e49 464
2fa65f52 465 if (TMath::Abs(particle->Eta()) < 1.5)
466 nMCTracks15 += particleWeight;
a9017e49 467
2fa65f52 468 if (TMath::Abs(particle->Eta()) < 2.0)
469 nMCTracks20 += particleWeight;
a9017e49 470
2fa65f52 471 nMCTracksAll += particleWeight;
a9017e49 472
2fa65f52 473 if (fParticleCorrection[0] || fParticleSpecies)
474 {
475 Int_t id = -1;
476 switch (TMath::Abs(particle->GetPdgCode()))
477 {
478 case 211: id = 0; break;
479 case 321: id = 1; break;
480 case 2212: id = 2; break;
481 default: id = 3; break;
482 }
a9017e49 483
2fa65f52 484 if (TMath::Abs(particle->Eta()) < etaRange)
485 nMCTracksSpecies[id]++;
486
487 if (fParticleCorrection[id])
488 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
489 }
490 } // end of mc particle
a9017e49 491
2fa65f52 492 fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, (Int_t) nMCTracks05, (Int_t) nMCTracks09, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
a9017e49 493
2fa65f52 494 if (eventTriggered && eventVertex)
a9017e49 495 {
2fa65f52 496 Int_t nESDTracks05 = 0;
497 Int_t nESDTracks09 = 0;
498 Int_t nESDTracks15 = 0;
499 Int_t nESDTracks20 = 0;
a9017e49 500
2fa65f52 501 // tracks per particle species, in |eta| < 2 (systematic study)
502 Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
503 for (Int_t i = 0; i<7; ++i)
504 nESDTracksSpecies[i] = 0;
a9017e49 505
2fa65f52 506 Bool_t* foundPrimaries = new Bool_t[nPrim]; // to prevent double counting
507 for (Int_t i=0; i<nPrim; i++)
508 foundPrimaries[i] = kFALSE;
a9017e49 509
2fa65f52 510 Bool_t* foundTracks = new Bool_t[nMCPart]; // to prevent double counting
511 for (Int_t i=0; i<nMCPart; i++)
512 foundTracks[i] = kFALSE;
a9017e49 513
2fa65f52 514 for (Int_t i=0; i<inputCount; ++i)
515 {
516 Float_t eta = etaArr[i];
517 Int_t label = labelArr[i];
a9017e49 518
2fa65f52 519 Int_t particleWeight = 1;
520
521 // in case of systematic study, weight according to the change of the pt spectrum
522 if (fPtSpectrum)
523 {
524 TParticle* mother = 0;
525
526 // preserve label for later
527 Int_t labelCopy = label;
528 if (labelCopy >= 0)
529 labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
530 if (labelCopy >= 0)
531 mother = stack->Particle(labelCopy);
532
533 // in case of pt study we do not count particles w/o label, because they cannot be scaled
534 if (!mother)
535 continue;
536
537 // it cannot be just multiplied because we cannot count "half of a particle"
538 // instead a random generator decides if the particle is counted twice (if value > 1)
539 // or not (if value < 0)
540 Int_t bin = fPtSpectrum->FindBin(mother->Pt());
541 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
542 {
543 Float_t factor = fPtSpectrum->GetBinContent(bin);
544 if (factor > 0)
545 {
546 Float_t random = gRandom->Uniform();
547 if (factor > 1 && random < factor - 1)
548 {
549 particleWeight = 2;
550 }
551 else if (factor < 1 && random < 1 - factor)
552 particleWeight = 0;
553 }
554 }
555 }
556
557 //Printf("ESD weight is: %d", particleWeight);
558
559 if (TMath::Abs(eta) < 0.5)
560 nESDTracks05 += particleWeight;
a9017e49 561
2fa65f52 562 if (TMath::Abs(eta) < 0.9)
563 nESDTracks09 += particleWeight;
a9017e49 564
2fa65f52 565 if (TMath::Abs(eta) < 1.5)
566 nESDTracks15 += particleWeight;
a9017e49 567
2fa65f52 568 if (TMath::Abs(eta) < 2.0)
569 nESDTracks20 += particleWeight;
a9017e49 570
a9017e49 571
2fa65f52 572 if (fParticleCorrection[0] || fParticleSpecies)
573 {
574 Int_t motherLabel = -1;
575 TParticle* mother = 0;
576
577 // find mother
578 if (label >= 0)
579 motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
580 if (motherLabel >= 0)
581 mother = stack->Particle(motherLabel);
582
583 if (!mother)
584 {
585 // count tracks that did not have a label
586 if (TMath::Abs(eta) < etaRange)
587 nESDTracksSpecies[4]++;
588 continue;
589 }
590
591 // get particle type (pion, proton, kaon, other)
592 Int_t id = -1;
593 switch (TMath::Abs(mother->GetPdgCode()))
594 {
595 case 211: id = 0; break;
596 case 321: id = 1; break;
597 case 2212: id = 2; break;
598 default: id = 3; break;
599 }
600
601 // double counting is ok for particle ratio study
602 if (TMath::Abs(eta) < etaRange)
603 nESDTracksSpecies[id]++;
604
605 // double counting is not ok for efficiency study
606
607 // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
608 if (foundTracks[label])
609 {
610 if (TMath::Abs(eta) < etaRange)
611 nESDTracksSpecies[6]++;
612 continue;
613 }
614 foundTracks[label] = kTRUE;
615
616 // particle (primary) already counted?
617 if (foundPrimaries[motherLabel])
618 {
619 if (TMath::Abs(eta) < etaRange)
620 nESDTracksSpecies[5]++;
621 continue;
622 }
623 foundPrimaries[motherLabel] = kTRUE;
624
625 if (fParticleCorrection[id])
626 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
627 }
628 }
629
630 delete[] foundTracks;
631 delete[] foundPrimaries;
a9017e49 632
2fa65f52 633 if ((Int_t) nMCTracks20 == 0 && nESDTracks20 > 3)
634 printf("WARNING: Event has %d generated and %d reconstructed...\n", nMCTracks20, nESDTracks20);
a9017e49 635
2fa65f52 636 // fill response matrix using vtxMC (best guess)
637 fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks09, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
a9017e49 638
2fa65f52 639 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
640
641 if (fParticleSpecies)
642 fParticleSpecies->Fill(vtxMC[2], nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
643 }
644 }
645 }
a9017e49 646
2fa65f52 647 delete etaArr;
648 delete labelArr;
a9017e49 649}
650
651void AliMultiplicityTask::Terminate(Option_t *)
652{
653 // The Terminate() function is the last function to be called during
654 // a query. It always runs on the client, it can be used to present
655 // the results graphically or save the results to file.
656
657 fOutput = dynamic_cast<TList*> (GetOutputData(0));
658 if (!fOutput) {
659 Printf("ERROR: fOutput not available");
660 return;
661 }
662
663 fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
664 for (Int_t i = 0; i < 4; ++i)
665 fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
666 fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
667
668 if (!fMultiplicity)
669 {
670 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
671 return;
672 }
673
674 TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
675
676 fMultiplicity->SaveHistograms();
677 for (Int_t i = 0; i < 4; ++i)
678 if (fParticleCorrection[i])
679 fParticleCorrection[i]->SaveHistograms();
680 if (fParticleSpecies)
681 fParticleSpecies->Write();
682
683 TObjString option(fOption);
684 option.Write();
685
686 file->Close();
687
688 Printf("Writting result to multiplicityMC.root");
689}