update to comply with new GetLabel function
[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
31#include "esdTrackCuts/AliESDtrackCuts.h"
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),
2fa65f52 43 fAnalysisMode(kSPD),
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
2fa65f52 98 if (fAnalysisMode == kSPD)
a9017e49 99 tree->SetBranchStatus("fSPDMult*", 1);
a9017e49 100
2fa65f52 101 if (fAnalysisMode == kTPC) {
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 {
200 AliDebug(AliLog::kError, "ESD branch not available");
201 return;
202 }
203
2fa65f52 204 // MB1 definition
205 //Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1);
206 // only FASTOR
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;
225 if (fAnalysisMode == 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 }
253 else if (fAnalysisMode == kTPC)
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) {
289 case kTPC: etaRange = 0.9; break;
290 case kSPD: etaRange = 2.0; break;
291 }
a9017e49 292
2fa65f52 293 if (!fReadMC) // Processing of ESD information
294 {
295 if (eventTriggered && eventVertex)
296 {
297 Int_t nESDTracks05 = 0;
298 Int_t nESDTracks09 = 0;
299 Int_t nESDTracks15 = 0;
300 Int_t nESDTracks20 = 0;
a9017e49 301
2fa65f52 302 for (Int_t i=0; i<inputCount; ++i)
303 {
304 Float_t eta = etaArr[i];
a9017e49 305
2fa65f52 306 if (TMath::Abs(eta) < 0.5)
307 nESDTracks05++;
a9017e49 308
2fa65f52 309 if (TMath::Abs(eta) < 0.9)
310 nESDTracks09++;
a9017e49 311
2fa65f52 312 if (TMath::Abs(eta) < 1.5)
313 nESDTracks15++;
314
315 if (TMath::Abs(eta) < 2.0)
316 nESDTracks20++;
a9017e49 317 }
2fa65f52 318
319 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
320 }
321 }
322 else if (fReadMC) // Processing of MC information
323 {
324 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
325 if (!eventHandler) {
326 Printf("ERROR: Could not retrieve MC event handler");
327 return;
a9017e49 328 }
a9017e49 329
2fa65f52 330 AliMCEvent* mcEvent = eventHandler->MCEvent();
331 if (!mcEvent) {
332 Printf("ERROR: Could not retrieve MC event");
333 return;
334 }
a9017e49 335
2fa65f52 336 AliStack* stack = mcEvent->Stack();
337 if (!stack)
338 {
339 AliDebug(AliLog::kError, "Stack not available");
340 return;
341 }
342
343 AliHeader* header = mcEvent->Header();
344 if (!header)
345 {
346 AliDebug(AliLog::kError, "Header not available");
347 return;
348 }
a9017e49 349
2fa65f52 350 Bool_t processEvent = kTRUE;
351 if (fSelectProcessType > 0)
352 {
353 // getting process information; NB: this only works for Pythia
354 Int_t processtype = AliPWG0Helper::GetPythiaEventProcessType(header);
a9017e49 355
2fa65f52 356 processEvent = kFALSE;
a9017e49 357
2fa65f52 358 // non diffractive
359 if (fSelectProcessType == 1 && processtype !=92 && processtype !=93 && processtype != 94)
360 processEvent = kTRUE;
a9017e49 361
2fa65f52 362 // single diffractive
363 if (fSelectProcessType == 2 && (processtype == 92 || processtype == 93))
364 processEvent = kTRUE;
a9017e49 365
2fa65f52 366 // double diffractive
367 if (fSelectProcessType == 3 && processtype == 94)
368 processEvent = kTRUE;
a9017e49 369
2fa65f52 370 if (!processEvent)
371 AliDebug(AliLog::kDebug, Form("Skipping this event, because it is not of the requested process type (%d)", processtype));
372 }
a9017e49 373
374 // systematic study: 10% lower efficiency
375 if (fSystSkipParticles && (gRandom->Uniform() < 0.1))
2fa65f52 376 processEvent = kFALSE;
a9017e49 377
2fa65f52 378 if (processEvent)
379 {
380 // get the MC vertex
381 AliGenEventHeader* genHeader = header->GenEventHeader();
382 TArrayF vtxMC(3);
383 genHeader->PrimaryVertex(vtxMC);
384
385 // get number of tracks from MC
386 // loop over mc particles
387 Int_t nPrim = stack->GetNprimary();
388 Int_t nMCPart = stack->GetNtrack();
389
390 // tracks in different eta ranges
391 Int_t nMCTracks05 = 0;
392 Int_t nMCTracks09 = 0;
393 Int_t nMCTracks15 = 0;
394 Int_t nMCTracks20 = 0;
395 Int_t nMCTracksAll = 0;
396
397 // tracks per particle species, in |eta| < 2 (systematic study)
398 Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
399 for (Int_t i = 0; i<4; ++i)
400 nMCTracksSpecies[i] = 0;
401
402 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
403 {
404 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
a9017e49 405
2fa65f52 406 TParticle* particle = stack->Particle(iMc);
a9017e49 407
2fa65f52 408 if (!particle)
409 {
410 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
411 continue;
412 }
a9017e49 413
2fa65f52 414 Bool_t debug = kFALSE;
415 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
416 {
417 //printf("%d) DROPPED ", iMC);
418 //particle->Print();
419 continue;
420 }
a9017e49 421
2fa65f52 422 //printf("%d) OK ", iMC);
423 //particle->Print();
a9017e49 424
2fa65f52 425 //if (particle->Pt() < kPtCut)
426 // continue;
a9017e49 427
2fa65f52 428 Int_t particleWeight = 1;
a9017e49 429
2fa65f52 430 Float_t pt = particle->Pt();
a9017e49 431
2fa65f52 432 // in case of systematic study, weight according to the change of the pt spectrum
433 // it cannot be just multiplied because we cannot count "half of a particle"
434 // instead a random generator decides if the particle is counted twice (if value > 1)
435 // or not (if value < 0)
436 if (fPtSpectrum)
a9017e49 437 {
2fa65f52 438 Int_t bin = fPtSpectrum->FindBin(pt);
439 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
a9017e49 440 {
2fa65f52 441 Float_t factor = fPtSpectrum->GetBinContent(bin);
442 if (factor > 0)
443 {
444 Float_t random = gRandom->Uniform();
445 if (factor > 1 && random < factor - 1)
446 {
447 particleWeight = 2;
448 }
449 else if (factor < 1 && random < 1 - factor)
450 particleWeight = 0;
451 }
a9017e49 452 }
a9017e49 453 }
a9017e49 454
2fa65f52 455 //Printf("MC weight is: %d", particleWeight);
a9017e49 456
2fa65f52 457 if (TMath::Abs(particle->Eta()) < 0.5)
458 nMCTracks05 += particleWeight;
a9017e49 459
2fa65f52 460 if (TMath::Abs(particle->Eta()) < 0.9)
461 nMCTracks09 += particleWeight;
a9017e49 462
2fa65f52 463 if (TMath::Abs(particle->Eta()) < 1.5)
464 nMCTracks15 += particleWeight;
a9017e49 465
2fa65f52 466 if (TMath::Abs(particle->Eta()) < 2.0)
467 nMCTracks20 += particleWeight;
a9017e49 468
2fa65f52 469 nMCTracksAll += particleWeight;
a9017e49 470
2fa65f52 471 if (fParticleCorrection[0] || fParticleSpecies)
472 {
473 Int_t id = -1;
474 switch (TMath::Abs(particle->GetPdgCode()))
475 {
476 case 211: id = 0; break;
477 case 321: id = 1; break;
478 case 2212: id = 2; break;
479 default: id = 3; break;
480 }
a9017e49 481
2fa65f52 482 if (TMath::Abs(particle->Eta()) < etaRange)
483 nMCTracksSpecies[id]++;
484
485 if (fParticleCorrection[id])
486 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
487 }
488 } // end of mc particle
a9017e49 489
2fa65f52 490 fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, (Int_t) nMCTracks05, (Int_t) nMCTracks09, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
a9017e49 491
2fa65f52 492 if (eventTriggered && eventVertex)
a9017e49 493 {
2fa65f52 494 Int_t nESDTracks05 = 0;
495 Int_t nESDTracks09 = 0;
496 Int_t nESDTracks15 = 0;
497 Int_t nESDTracks20 = 0;
a9017e49 498
2fa65f52 499 // tracks per particle species, in |eta| < 2 (systematic study)
500 Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
501 for (Int_t i = 0; i<7; ++i)
502 nESDTracksSpecies[i] = 0;
a9017e49 503
2fa65f52 504 Bool_t* foundPrimaries = new Bool_t[nPrim]; // to prevent double counting
505 for (Int_t i=0; i<nPrim; i++)
506 foundPrimaries[i] = kFALSE;
a9017e49 507
2fa65f52 508 Bool_t* foundTracks = new Bool_t[nMCPart]; // to prevent double counting
509 for (Int_t i=0; i<nMCPart; i++)
510 foundTracks[i] = kFALSE;
a9017e49 511
2fa65f52 512 for (Int_t i=0; i<inputCount; ++i)
513 {
514 Float_t eta = etaArr[i];
515 Int_t label = labelArr[i];
a9017e49 516
2fa65f52 517 Int_t particleWeight = 1;
518
519 // in case of systematic study, weight according to the change of the pt spectrum
520 if (fPtSpectrum)
521 {
522 TParticle* mother = 0;
523
524 // preserve label for later
525 Int_t labelCopy = label;
526 if (labelCopy >= 0)
527 labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
528 if (labelCopy >= 0)
529 mother = stack->Particle(labelCopy);
530
531 // in case of pt study we do not count particles w/o label, because they cannot be scaled
532 if (!mother)
533 continue;
534
535 // it cannot be just multiplied because we cannot count "half of a particle"
536 // instead a random generator decides if the particle is counted twice (if value > 1)
537 // or not (if value < 0)
538 Int_t bin = fPtSpectrum->FindBin(mother->Pt());
539 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
540 {
541 Float_t factor = fPtSpectrum->GetBinContent(bin);
542 if (factor > 0)
543 {
544 Float_t random = gRandom->Uniform();
545 if (factor > 1 && random < factor - 1)
546 {
547 particleWeight = 2;
548 }
549 else if (factor < 1 && random < 1 - factor)
550 particleWeight = 0;
551 }
552 }
553 }
554
555 //Printf("ESD weight is: %d", particleWeight);
556
557 if (TMath::Abs(eta) < 0.5)
558 nESDTracks05 += particleWeight;
a9017e49 559
2fa65f52 560 if (TMath::Abs(eta) < 0.9)
561 nESDTracks09 += particleWeight;
a9017e49 562
2fa65f52 563 if (TMath::Abs(eta) < 1.5)
564 nESDTracks15 += particleWeight;
a9017e49 565
2fa65f52 566 if (TMath::Abs(eta) < 2.0)
567 nESDTracks20 += particleWeight;
a9017e49 568
a9017e49 569
2fa65f52 570 if (fParticleCorrection[0] || fParticleSpecies)
571 {
572 Int_t motherLabel = -1;
573 TParticle* mother = 0;
574
575 // find mother
576 if (label >= 0)
577 motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
578 if (motherLabel >= 0)
579 mother = stack->Particle(motherLabel);
580
581 if (!mother)
582 {
583 // count tracks that did not have a label
584 if (TMath::Abs(eta) < etaRange)
585 nESDTracksSpecies[4]++;
586 continue;
587 }
588
589 // get particle type (pion, proton, kaon, other)
590 Int_t id = -1;
591 switch (TMath::Abs(mother->GetPdgCode()))
592 {
593 case 211: id = 0; break;
594 case 321: id = 1; break;
595 case 2212: id = 2; break;
596 default: id = 3; break;
597 }
598
599 // double counting is ok for particle ratio study
600 if (TMath::Abs(eta) < etaRange)
601 nESDTracksSpecies[id]++;
602
603 // double counting is not ok for efficiency study
604
605 // 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)
606 if (foundTracks[label])
607 {
608 if (TMath::Abs(eta) < etaRange)
609 nESDTracksSpecies[6]++;
610 continue;
611 }
612 foundTracks[label] = kTRUE;
613
614 // particle (primary) already counted?
615 if (foundPrimaries[motherLabel])
616 {
617 if (TMath::Abs(eta) < etaRange)
618 nESDTracksSpecies[5]++;
619 continue;
620 }
621 foundPrimaries[motherLabel] = kTRUE;
622
623 if (fParticleCorrection[id])
624 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
625 }
626 }
627
628 delete[] foundTracks;
629 delete[] foundPrimaries;
a9017e49 630
2fa65f52 631 if ((Int_t) nMCTracks20 == 0 && nESDTracks20 > 3)
632 printf("WARNING: Event has %d generated and %d reconstructed...\n", nMCTracks20, nESDTracks20);
a9017e49 633
2fa65f52 634 // fill response matrix using vtxMC (best guess)
635 fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks09, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
a9017e49 636
2fa65f52 637 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
638
639 if (fParticleSpecies)
640 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]);
641 }
642 }
643 }
a9017e49 644
2fa65f52 645 delete etaArr;
646 delete labelArr;
a9017e49 647}
648
649void AliMultiplicityTask::Terminate(Option_t *)
650{
651 // The Terminate() function is the last function to be called during
652 // a query. It always runs on the client, it can be used to present
653 // the results graphically or save the results to file.
654
655 fOutput = dynamic_cast<TList*> (GetOutputData(0));
656 if (!fOutput) {
657 Printf("ERROR: fOutput not available");
658 return;
659 }
660
661 fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
662 for (Int_t i = 0; i < 4; ++i)
663 fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
664 fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
665
666 if (!fMultiplicity)
667 {
668 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
669 return;
670 }
671
672 TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
673
674 fMultiplicity->SaveHistograms();
675 for (Int_t i = 0; i < 4; ++i)
676 if (fParticleCorrection[i])
677 fParticleCorrection[i]->SaveHistograms();
678 if (fParticleSpecies)
679 fParticleSpecies->Write();
680
681 TObjString option(fOption);
682 option.Write();
683
684 file->Close();
685
686 Printf("Writting result to multiplicityMC.root");
687}