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