2c413cea6ec4d5c51d586f72672d9145a46da9df
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / AliMultiplicityTask.cxx
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 <TH1D.h>
12 #include <TH2F.h>
13 #include <TH3F.h>
14 #include <TParticle.h>
15 #include <TRandom.h>
16 #include <TNtuple.h>
17 #include <TObjString.h>
18 #include <TF1.h>
19
20 #include <AliLog.h>
21 #include <AliESDVertex.h>
22 #include <AliESDEvent.h>
23 #include <AliStack.h>
24 #include <AliHeader.h>
25 #include <AliGenEventHeader.h>
26 #include <AliMultiplicity.h>
27 #include <AliAnalysisManager.h>
28 #include <AliMCEventHandler.h>
29 #include <AliMCEvent.h>
30 #include <AliESDInputHandler.h>
31
32 #include "AliESDtrackCuts.h"
33 #include "AliPWG0Helper.h"
34 #include "AliMultiplicityCorrection.h"
35 #include "AliCorrection.h"
36 #include "AliCorrectionMatrix3D.h"
37 #include "AliPhysicsSelection.h"
38 #include "AliTriggerAnalysis.h"
39 #include "AliVEvent.h"
40
41 ClassImp(AliMultiplicityTask)
42
43 AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
44   AliAnalysisTask("AliMultiplicityTask", ""),
45   fESD(0),
46   fOption(opt),
47   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn)),
48   fTrigger(AliTriggerAnalysis::kMB1),
49   fDeltaPhiCut(-1),
50   fDiffTreatment(AliPWG0Helper::kMCFlags),
51   fReadMC(kFALSE),
52   fUseMCVertex(kFALSE),
53   fMultiplicity(0),
54   fEsdTrackCuts(0),
55   fSystSkipParticles(kFALSE),
56   fSelectProcessType(0),
57   fParticleSpecies(0),
58   fdNdpT(0),
59   fPtSpectrum(0),
60   fTemp1(0),
61   fTemp2(0),
62   fVertex(0),
63   fEtaPhi(0),
64   fOutput(0)
65 {
66   //
67   // Constructor. Initialization of pointers
68   //
69
70   for (Int_t i = 0; i<8; i++)
71     fParticleCorrection[i] = 0;
72     
73   for (Int_t i=0; i<3; i++)
74     fEta[i] = 0;
75
76   // Define input and output slots here
77   DefineInput(0, TChain::Class());
78   DefineOutput(0, TList::Class());
79
80   if (fOption.Contains("only-process-type-nd"))
81     fSelectProcessType = 1;
82
83   if (fOption.Contains("only-process-type-sd"))
84     fSelectProcessType = 2;
85
86   if (fOption.Contains("only-process-type-dd"))
87     fSelectProcessType = 3;
88
89   if (fSelectProcessType != 0)
90     AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
91 }
92
93 AliMultiplicityTask::~AliMultiplicityTask()
94 {
95   //
96   // Destructor
97   //
98
99   // histograms are in the output list and deleted when the output
100   // list is deleted by the TSelector dtor
101
102   if (fOutput&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
103     delete fOutput;
104     fOutput = 0;
105   }
106 }
107
108 //________________________________________________________________________
109 void AliMultiplicityTask::ConnectInputData(Option_t *)
110 {
111   // Connect ESD
112   // Called once
113
114   Printf("AliMultiplicityTask::ConnectInputData called");
115
116   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
117   if (!tree) {
118     Printf("ERROR: Could not read tree from input slot 0");
119   } else {
120     // Disable all branches and enable only the needed ones
121     /*
122     tree->SetBranchStatus("*", 0);
123
124     tree->SetBranchStatus("AliESDHeader*", 1);
125     tree->SetBranchStatus("*Vertex*", 1);
126
127     if (fAnalysisMode & AliPWG0Helper::kSPD) {
128       tree->SetBranchStatus("AliMultiplicity*", 1);
129     }
130
131     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
132       //AliESDtrackCuts::EnableNeededBranches(tree);
133       tree->SetBranchStatus("Tracks*", 1);
134     }
135     */
136
137     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
138
139     if (!esdH) {
140       Printf("ERROR: Could not get ESDInputHandler");
141     } else
142       fESD = esdH->GetEvent();
143   }
144
145   // disable info messages of AliMCEvent (per event)
146   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
147 }
148
149 void AliMultiplicityTask::CreateOutputObjects()
150 {
151   // create result objects and add to output list
152
153   fOutput = new TList;
154   fOutput->SetOwner();
155
156   fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
157   fOutput->Add(fMultiplicity);
158   
159   fdNdpT = new TH1F("fdNdpT", "fdNdpT", 1000, 0, 10);
160   fdNdpT->Sumw2();
161   fOutput->Add(fdNdpT);
162
163   if (fOption.Contains("particle-efficiency"))
164     for (Int_t i = 0; i<8; i++)
165     {
166       fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
167       fOutput->Add(fParticleCorrection[i]);
168     }
169
170   if (fOption.Contains("pt-spectrum-hist"))
171   {
172     TFile* file = TFile::Open("ptspectrum_fit.root");
173     if (file)
174     {
175       TString subStr(fOption(fOption.Index("pt-spectrum")+17, 3));
176       TString histName(Form("ptspectrum_%s", subStr.Data()));
177       AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
178       fPtSpectrum = (TH1D*) file->Get(histName);
179       if (!fPtSpectrum)   
180         AliError("Histogram not found");
181     }
182     else
183       AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
184   }
185
186   if (fOption.Contains("pt-spectrum-func"))
187   {
188     if (fPtSpectrum)
189     {
190       Printf("Using function for systematic p_t study");
191     }
192     else
193     {
194       Printf("ERROR: Could not find function for systematic p_t study");
195       fPtSpectrum = new TH1D("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
196       fPtSpectrum->SetBinContent(1, 1);
197     }
198   }
199
200   if (fPtSpectrum)
201     Printf("WARNING: Systematic study enabled. Pt spectrum will be modified");
202   
203   if (fOption.Contains("particle-species")) {
204     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");
205     fOutput->Add(fParticleSpecies);
206   }
207
208   fTemp1 = new TH2F("fTemp1", "fTemp1", 100, -0.5, 99.5, 100, -0.5, 99.5);
209   fOutput->Add(fTemp1);
210   
211   for (Int_t i=0; i<3; i++)
212   {
213     fEta[i] = new TH1F(Form("fEta_%d", i), ";#eta", 100, -2, 2);
214     fOutput->Add(fEta[i]);
215   }
216   
217   fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
218   fOutput->Add(fVertex);
219   
220   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
221   fOutput->Add(fEtaPhi);
222   
223   // TODO set seed for random generator
224 }
225
226 void AliMultiplicityTask::Exec(Option_t*)
227 {
228   // process the event
229
230   // Check prerequisites
231   if (!fESD)
232   {
233     AliDebug(AliLog::kError, "ESD not available");
234     return;
235   }
236
237   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
238   if (!inputHandler)
239   {
240     Printf("ERROR: Could not receive input handler");
241     return;
242   }
243     
244   Bool_t eventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
245
246   static AliTriggerAnalysis* triggerAnalysis = 0;
247   if (!triggerAnalysis)
248   {
249     AliPhysicsSelection* physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
250     triggerAnalysis = physicsSelection->GetTriggerAnalysis();
251   }
252   if (eventTriggered)
253     eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
254     
255   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
256   if (vtxESD && !AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
257     vtxESD = 0;
258     
259   // remove vertices outside +- 15 cm
260   if (vtxESD && TMath::Abs(vtxESD->GetZv()) > 15)
261     vtxESD = 0;
262   
263   Bool_t eventVertex = (vtxESD != 0);
264
265   Double_t vtx[3];
266   if (vtxESD)
267   {
268     vtxESD->GetXYZ(vtx);
269     fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
270   }
271   
272   // post the data already here
273   PostData(0, fOutput);
274   
275   //const Float_t kPtCut = 0.3;
276
277   // create list of (label, eta) tuples
278   Int_t inputCount = 0;
279   Int_t* labelArr = 0;
280   Float_t* etaArr = 0;
281   Int_t nGoodTracks = 0;  // number of total good tracks is needed both in SPD and TPC loops if the mode is kTPCSPD
282   Int_t nESDTracks  = 0; // Total number of ESD tracks (including those not passing the cuts);
283   const int kRejBit = BIT(15); // set this bit in ESD tracks if it is rejected by a cut
284
285   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS || fAnalysisMode & AliPWG0Helper::kTPCSPD)
286   {
287     if (!fEsdTrackCuts)
288     {
289       AliDebug(AliLog::kError, "fESDTrackCuts not available");
290       return;
291     }
292
293     if (vtxESD)
294     {
295       // get multiplicity from ESD tracks
296       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
297       nGoodTracks = list->GetEntries();
298   
299       Int_t arraySize = nGoodTracks; // if we're also using tracklets, we need to increase this
300       if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
301         // if I have to also count clusters later on, I have to make sure the arrays are big enough
302         // Moreover, we need to keep track of rejected tracks to avoid double-counting
303         Printf( "Processing tracks");
304         const AliMultiplicity* mult = fESD->GetMultiplicity();
305         if (!mult)
306           {
307             AliDebug(AliLog::kError, "AliMultiplicity not available");
308             return;
309           }
310         arraySize += mult->GetNumberOfTracklets();
311         // allocate and fill array for rejected tracks
312         nESDTracks = fESD->GetNumberOfTracks();
313         for(Int_t itracks = 0; itracks < nESDTracks; itracks++){
314           AliESDtrack* track = fESD->GetTrack(itracks);
315           if( itracks!=track->GetID() ){
316             AliFatal("Track ID not matching track index!");
317           }
318           if (!fEsdTrackCuts->AcceptTrack(track)) {         
319             track->SetBit(kRejBit);
320           }
321         }
322       }
323
324       labelArr = new Int_t[arraySize];
325       etaArr = new Float_t[arraySize];
326   
327       // loop over esd tracks
328       for (Int_t i=0; i<nGoodTracks; i++)
329       {
330         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
331         if (!esdTrack)
332         {
333           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
334           continue;
335         }
336         
337         if (esdTrack->Pt() < 0.15){
338           esdTrack->SetBit(kRejBit);
339           continue;
340         }
341
342         if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
343           // Cuts by ruben to flag secondaries
344           if (esdTrack->IsOn(AliESDtrack::kMultSec) ) continue;
345           if (esdTrack->IsOn(AliESDtrack::kMultInV0)) continue;
346         }
347         
348         Float_t d0z0[2],covd0z0[3];
349         esdTrack->GetImpactParameters(d0z0,covd0z0);
350         Float_t sigma= 0.0050+0.0060/TMath::Power(esdTrack->Pt(),0.9);
351         Float_t d0max = 7.*sigma;
352         if (TMath::Abs(d0z0[0]) > d0max) {
353           //      rejLabelArr[irejCount++] = esdTrack->GetID(); // We
354           //      don't count the tracklet if the track is a
355           //      secondary, so this must be commented out
356           continue;
357         }
358   
359         if (vtxESD && TMath::Abs(vtx[2]) < 10)
360           fEtaPhi->Fill(esdTrack->Eta(), esdTrack->Phi());
361         
362         etaArr[inputCount] = esdTrack->Eta();
363         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
364         ++inputCount;
365       }
366   
367       delete list;
368     }
369   }
370   // SPD tracklets are used either in SPD standalone mode or in TPC+SPD mode (tracks + tracklets not using any cluster already used in tracks)
371   if (fAnalysisMode & AliPWG0Helper::kSPD || fAnalysisMode & AliPWG0Helper::kTPCSPD)
372   {
373
374     if (vtxESD)
375     {
376       AliDebug(AliLog::kInfo, "Processing tracklets");
377       // get tracklets
378       const AliMultiplicity* mult = fESD->GetMultiplicity();
379       if (!mult)
380       {
381         AliDebug(AliLog::kError, "AliMultiplicity not available");
382         if (labelArr) delete[] labelArr;
383         if (etaArr) delete[] etaArr;
384         return;
385       }
386   
387       Bool_t foundInEta10 = kFALSE;
388       if ((fAnalysisMode & AliPWG0Helper::kSPD) && !(fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS || fAnalysisMode & AliPWG0Helper::kTPCSPD)) {
389         // if we are counting both tracks and tracklets, these arrays were already initialized above 
390         AliDebug(AliLog::kInfo, "Booking arrays");
391         if (!labelArr) labelArr = new Int_t[mult->GetNumberOfTracklets()];
392         if (!etaArr) etaArr = new Float_t[mult->GetNumberOfTracklets()];
393       }
394       if (inputCount) foundInEta10 = kTRUE; // by definition, if we found a track.
395       
396       // get multiplicity from ITS tracklets
397       for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
398       {
399         //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
400   
401         Float_t deltaPhi = mult->GetDeltaPhi(i);
402         
403         if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
404           continue;
405   
406         if (fSystSkipParticles && gRandom->Uniform() < (0.0153))
407           {
408           Printf("Skipped tracklet!");
409           continue;
410           }
411         
412         // if counting tracks+tracklets, check if clusters where already used in tracks
413         if (fAnalysisMode & AliPWG0Helper::kTPCSPD) { 
414             Int_t id1,id2;
415             mult->GetTrackletTrackIDs(i,0,id1,id2);
416             if ( (id1>=0&& !fESD->GetTrack(id1)->TestBit(kRejBit) ) ||
417                  (id2>=0&& !fESD->GetTrack(id2)->TestBit(kRejBit) )
418                  ) {
419               printf ("Skipping tracklet: already used in track");  
420               continue; 
421             }
422         // Ruben also had this, but we're not using pure ITSSA tracks here:
423         // if (fSPDMult->FreeClustersTracklet(itr,1)) trITSSApure++; // not used in ITS_SA_Pure track
424         }
425                
426
427         etaArr[inputCount] = mult->GetEta(i);
428         if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
429         {
430           labelArr[inputCount] = mult->GetLabel(i, 0);
431         }
432         else
433           labelArr[inputCount] = -1;
434           
435         for (Int_t j=0; j<3; j++)
436         {
437           if (vtx[2] > fMultiplicity->GetVertexBegin(j) && vtx[2] < fMultiplicity->GetVertexEnd(j))
438             fEta[j]->Fill(etaArr[inputCount]);
439         }
440         
441         if (vtxESD && TMath::Abs(vtx[2]) < 10)
442           fEtaPhi->Fill(etaArr[inputCount], mult->GetPhi(i));
443         
444           // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
445         if (TMath::Abs(etaArr[inputCount]) < 1)
446           foundInEta10 = kTRUE;
447           
448         ++inputCount;
449       }
450       
451       if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
452         eventTriggered = kFALSE;
453     }
454   }
455
456   // eta range for nMCTracksSpecies, nESDTracksSpecies
457   Float_t etaRange = 1.0;
458 //   switch (fAnalysisMode) {
459 //     case AliPWG0Helper::kInvalid: break;
460 //     case AliPWG0Helper::kTPC:
461 //     case AliPWG0Helper::kTPCITS:
462 //      etaRange = 1.0; break;
463 //     case AliPWG0Helper::kSPD: etaRange = 1.0; break;
464 //     default: break;
465 //   }
466
467   if (!fReadMC) // Processing of ESD information
468   {
469     Int_t nESDTracks05 = 0;
470     Int_t nESDTracks10 = 0;
471     Int_t nESDTracks14 = 0;
472     
473     for (Int_t i=0; i<inputCount; ++i)
474     {
475       Float_t eta = etaArr[i];
476
477       if (TMath::Abs(eta) < 0.5)
478         nESDTracks05++;
479
480       if (TMath::Abs(eta) < 1.0)
481         nESDTracks10++;
482
483       if (TMath::Abs(eta) < 1.3)
484         nESDTracks14++;
485     }
486     
487     //if (nESDTracks05 >= 20 || nESDTracks10 >= 30 || nESDTracks14 >= 32)
488     //  Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d; Tracks: %d %d %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber(), nESDTracks05, nESDTracks10, nESDTracks14);
489
490     if (eventTriggered)
491       fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
492     
493     if (eventTriggered && eventVertex)
494       fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
495   }
496   else if (fReadMC)   // Processing of MC information
497   {
498     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
499     if (!eventHandler) {
500       Printf("ERROR: Could not retrieve MC event handler");
501       if (labelArr) delete[] labelArr;
502       if (etaArr) delete[] etaArr;
503       return;
504     }
505
506     AliMCEvent* mcEvent = eventHandler->MCEvent();
507     if (!mcEvent) {
508       Printf("ERROR: Could not retrieve MC event");
509       if (labelArr) delete[] labelArr;
510       if (etaArr) delete[] etaArr;
511       return;
512     }
513
514     AliStack* stack = mcEvent->Stack();
515     if (!stack)
516     {
517       AliDebug(AliLog::kError, "Stack not available");
518       if (labelArr) delete[] labelArr;
519       if (etaArr) delete[] etaArr;
520       return;
521     }
522     
523     AliHeader* header = mcEvent->Header();
524     if (!header)
525     {
526       AliDebug(AliLog::kError, "Header not available");
527       if (labelArr) delete[] labelArr;
528       if (etaArr) delete[] etaArr;
529       return;
530     }
531
532     if (fUseMCVertex)
533     {
534       Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
535       // get the MC vertex
536       AliGenEventHeader* genHeader = header->GenEventHeader();
537       TArrayF vtxMC(3);
538       genHeader->PrimaryVertex(vtxMC);
539
540       vtx[2] = vtxMC[2];
541     }
542     
543     // get process information
544     AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
545
546     Bool_t processEvent = kTRUE;
547     if (fSelectProcessType > 0)
548     {
549       processEvent = kFALSE;
550
551       // non diffractive
552       if (fSelectProcessType == 1 && processType == AliPWG0Helper::kND)
553         processEvent = kTRUE;
554
555       // single diffractive
556       if (fSelectProcessType == 2 && processType == AliPWG0Helper::kSD)
557         processEvent = kTRUE;
558
559       // double diffractive
560       if (fSelectProcessType == 3 && processType == AliPWG0Helper::kDD)
561         processEvent = kTRUE;
562
563       if (!processEvent)
564         Printf("Skipping this event, because it is not of the requested process type (%d)", (Int_t) processType);
565     }
566
567     if (processEvent)
568     {
569       // get the MC vertex
570       AliGenEventHeader* genHeader = header->GenEventHeader();
571       TArrayF vtxMC(3);
572       genHeader->PrimaryVertex(vtxMC);
573
574       // get number of tracks from MC
575       // loop over mc particles
576       Int_t nPrim  = stack->GetNprimary();
577       Int_t nMCPart = stack->GetNtrack();
578
579       // tracks in different eta ranges
580       Int_t nMCTracks05 = 0;
581       Int_t nMCTracks10 = 0;
582       Int_t nMCTracks14 = 0;
583       Int_t nMCTracksAll = 0;
584
585       // tracks per particle species, in |eta| < 2 (systematic study)
586       Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
587       for (Int_t i = 0; i<4; ++i)
588         nMCTracksSpecies[i] = 0;
589
590       for (Int_t iMc = 0; iMc < nPrim; ++iMc)
591       {
592         AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
593
594         TParticle* particle = stack->Particle(iMc);
595
596         if (!particle)
597         {
598           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
599           continue;
600         }
601
602         Bool_t debug = kFALSE;
603         if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
604         {
605           //printf("%d) DROPPED ", iMC);
606           //particle->Print();
607           continue;
608         }
609
610         //printf("%d) OK      ", iMC);
611         //particle->Print();
612
613         //if (particle->Pt() < kPtCut)
614         //  continue;
615
616         Int_t particleWeight = 1;
617
618         Float_t pt = particle->Pt();
619
620         // in case of systematic study, weight according to the change of the pt spectrum
621         // it cannot be just multiplied because we cannot count "half of a particle"
622         // instead a random generator decides if the particle is counted twice (if value > 1)
623         // or not (if value < 0)
624         if (fPtSpectrum)
625         {
626           Int_t bin = fPtSpectrum->FindBin(pt);
627           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
628           {
629             Float_t factor = fPtSpectrum->GetBinContent(bin);
630             if (factor > 0)
631             {
632               Float_t random = gRandom->Uniform();
633               if (factor > 1 && random < factor - 1)
634               {
635                 particleWeight = 2;
636               }
637               else if (factor < 1 && random < 1 - factor)
638                 particleWeight = 0;
639             }
640           }
641         }
642
643         //Printf("MC weight is: %d", particleWeight);
644
645         if (TMath::Abs(particle->Eta()) < 0.5)
646           nMCTracks05 += particleWeight;
647
648         if (TMath::Abs(particle->Eta()) < 1.0)
649           nMCTracks10 += particleWeight;
650
651         if (TMath::Abs(particle->Eta()) < 1.3)
652           nMCTracks14 += particleWeight;
653
654         nMCTracksAll += particleWeight;
655         
656         if (particle->Pt() > 0 && TMath::Abs(particle->Eta()) < 1.0)
657           fdNdpT->Fill(particle->Pt(), 1.0 / particle->Pt());
658
659         if (fParticleCorrection[0] || fParticleSpecies)
660         {
661           Int_t id = -1;
662           switch (TMath::Abs(particle->GetPdgCode()))
663           {
664             case 211: id = 0; break;
665             case 321: id = 1; break;
666             case 2212: id = 2; break;
667             default: id = 3; break;
668           }
669
670           if (TMath::Abs(particle->Eta()) < etaRange)
671             nMCTracksSpecies[id]++;
672             
673           if (fParticleCorrection[id])
674             fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
675         }
676       } // end of mc particle
677
678       fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
679
680       // ESD processing
681       Int_t nESDTracks05 = 0;
682       Int_t nESDTracks10 = 0;
683       Int_t nESDTracks14 = 0;
684
685       // tracks per particle species, in |eta| < 2 (systematic study)
686       Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
687       for (Int_t i = 0; i<7; ++i)
688         nESDTracksSpecies[i] = 0;
689
690       Bool_t* foundPrimaries = new Bool_t[nMCPart];   // to prevent double counting
691       for (Int_t i=0; i<nPrim; i++)
692         foundPrimaries[i] = kFALSE;
693
694       Bool_t* foundPrimaries2 = new Bool_t[nMCPart];   // to prevent double counting
695       for (Int_t i=0; i<nPrim; i++)
696         foundPrimaries2[i] = kFALSE;
697
698       Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
699       for (Int_t i=0; i<nMCPart; i++)
700         foundTracks[i] = kFALSE;
701
702       for (Int_t i=0; i<inputCount; ++i)
703       {
704         Float_t eta = etaArr[i];
705         Int_t label = labelArr[i];
706
707         Int_t particleWeight = 1;
708
709         // in case of systematic study, weight according to the change of the pt spectrum
710         if (fPtSpectrum)
711         {
712           TParticle* mother = 0;
713
714           // preserve label for later
715           Int_t labelCopy = label;
716           if (labelCopy >= 0)
717             labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
718           if (labelCopy >= 0)
719             mother = stack->Particle(labelCopy);
720
721           // in case of pt study we do not count particles w/o label, because they cannot be scaled
722           if (!mother)
723             continue;
724
725           // it cannot be just multiplied because we cannot count "half of a particle"
726           // instead a random generator decides if the particle is counted twice (if value > 1) 
727           // or not (if value < 0)
728           Int_t bin = fPtSpectrum->FindBin(mother->Pt());
729           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
730           {
731             Float_t factor = fPtSpectrum->GetBinContent(bin);
732             if (factor > 0)
733             {
734               Float_t random = gRandom->Uniform();
735               if (factor > 1 && random < factor - 1)
736               {
737                 particleWeight = 2;
738               }
739               else if (factor < 1 && random < 1 - factor)
740                 particleWeight = 0;
741             }
742           }
743         }
744
745         //Printf("ESD weight is: %d", particleWeight);
746
747         if (TMath::Abs(eta) < 0.5)
748           nESDTracks05 += particleWeight;
749
750         if (TMath::Abs(eta) < 1.0)
751           nESDTracks10 += particleWeight;
752
753         if (TMath::Abs(eta) < 1.3)
754           nESDTracks14 += particleWeight;
755
756         if (fParticleSpecies)
757         {
758           Int_t motherLabel = -1;
759           TParticle* mother = 0;
760
761           // find mother
762           if (label >= 0)
763             motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
764           if (motherLabel >= 0)
765             mother = stack->Particle(motherLabel);
766
767           if (!mother)
768           {
769             // count tracks that did not have a label
770             if (TMath::Abs(eta) < etaRange)
771               nESDTracksSpecies[4]++;
772           }
773           else
774           {
775             // get particle type (pion, proton, kaon, other) of mother
776             Int_t idMother = -1;
777             switch (TMath::Abs(mother->GetPdgCode()))
778             {
779               case 211: idMother = 0; break;
780               case 321: idMother = 1; break;
781               case 2212: idMother = 2; break;
782               default: idMother = 3; break;
783             }
784
785             // double counting is ok for particle ratio study
786             if (TMath::Abs(eta) < etaRange)
787               nESDTracksSpecies[idMother]++;
788
789             // double counting is not ok for efficiency study
790
791             // 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)
792             if (foundTracks[label])
793             {
794               if (TMath::Abs(eta) < etaRange)
795                 nESDTracksSpecies[6]++;
796             }
797             else
798             {
799               foundTracks[label] = kTRUE;
800
801               // particle (primary) already counted?
802               if (foundPrimaries[motherLabel])
803               {
804                 if (TMath::Abs(eta) < etaRange)
805                   nESDTracksSpecies[5]++;
806               }
807               else
808                 foundPrimaries[motherLabel] = kTRUE;
809             }
810           }
811         }
812
813         if (fParticleCorrection[0])
814         {
815           if (label >= 0 && stack->IsPhysicalPrimary(label))
816           {
817             TParticle* particle = stack->Particle(label);
818
819             // get particle type (pion, proton, kaon, other)
820             Int_t id = -1;
821             switch (TMath::Abs(particle->GetPdgCode()))
822             {
823               case 211: id = 0; break;
824               case 321: id = 1; break;
825               case 2212: id = 2; break;
826               default: id = 3; break;
827             }
828
829             // todo check if values are not completely off??
830
831             // particle (primary) already counted?
832             if (!foundPrimaries2[label])
833             {
834               foundPrimaries2[label] = kTRUE;
835               fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
836             }
837           }
838         }
839       }
840         
841       if (fParticleCorrection[0])
842       {
843         // if the particle decays/stops before this radius we do not see it
844         // 8cm larger than SPD layer 2
845         // 123cm TPC radius where a track has about 50 clusters (cut limit)          
846         const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
847                 
848         // loop over all primaries that have not been found
849         for (Int_t i=0; i<nPrim; i++)
850         {
851           // already found
852           if (foundPrimaries2[i])
853             continue;
854             
855           TParticle* particle = 0;
856           TClonesArray* trackrefs = 0;
857           mcEvent->GetParticleAndTR(i, particle, trackrefs);
858           
859           // true primary and charged
860           if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
861             continue;              
862           
863           //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
864           if (TMath::Abs(particle->Eta()) > 3)
865             continue;
866           
867           // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
868           
869           // get particle type (pion, proton, kaon, other)
870           Int_t id = -1;
871           switch (TMath::Abs(particle->GetPdgCode()))
872           {
873             case 211: id = 4; break;
874             case 321: id = 5; break;
875             case 2212: id = 6; break;
876             default: id = 7; break;
877           }            
878           
879           if (!fParticleCorrection[id])
880             continue;
881             
882           // get last track reference
883           AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
884           
885           if (!trackref)
886           {
887             Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
888             particle->Print();
889             continue;
890           }
891             
892           // particle in tracking volume long enough...
893           if (trackref->R() > endRadius)
894             continue;  
895           
896           if (particle->GetLastDaughter() >= 0)
897           {
898             Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
899             //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
900             if (uID == kPDecay)
901             {
902               // decayed
903               
904               Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
905               particle->Print();
906               Printf("Daugthers:");
907               for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
908                 stack->Particle(d)->Print();
909               Printf(" ");
910               
911               fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
912               continue;
913             }
914           }
915           
916           if (trackref->DetectorId() == -1)
917           {
918             // stopped
919             Printf("Particle %d stopped at %f:", i, trackref->R());
920             particle->Print();
921             Printf(" ");
922             
923             fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
924             continue;
925           }
926           
927           Printf("Particle %d simply not tracked", i);
928           particle->Print();
929           Printf(" ");
930         }
931       }
932       
933       delete[] foundTracks;
934       delete[] foundPrimaries;
935       delete[] foundPrimaries2;
936
937 //         if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
938 //         {
939 //             TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
940 //             printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
941 //         }
942
943       if (eventTriggered)
944       {
945         fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
946         fMultiplicity->FillNoVertexEvent(vtxMC[2], eventVertex, nMCTracks05, nMCTracks10, nMCTracks14, nESDTracks05, nESDTracks10, nESDTracks14);
947 //         if (!eventVertex)
948 //         {
949 //           if (nESDTracks05 == 0)
950 //             fTemp1->Fill(nMCTracks05, nESDTracks05);
951 //         }
952       }
953       
954       if (eventTriggered && eventVertex)
955       {
956         // fill response matrix using vtxMC (best guess)
957         fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks14,  nMCTracksAll,  nESDTracks05, nESDTracks10, nESDTracks14);
958
959         fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
960
961         if (fParticleSpecies)
962           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]);
963       }
964     }
965   }
966
967   if (etaArr)
968     delete[] etaArr;
969   if (labelArr)
970     delete[] labelArr;
971
972 }
973
974 void AliMultiplicityTask::Terminate(Option_t *)
975 {
976   // The Terminate() function is the last function to be called during
977   // a query. It always runs on the client, it can be used to present
978   // the results graphically or save the results to file.
979
980   fOutput = dynamic_cast<TList*> (GetOutputData(0));
981   if (!fOutput) {
982     Printf("ERROR: fOutput not available");
983     return;
984   }
985
986   fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
987   for (Int_t i = 0; i < 8; ++i)
988     fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
989   fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
990   
991   fdNdpT = dynamic_cast<TH1*> (fOutput->FindObject("fdNdpT"));
992
993   if (!fMultiplicity)
994   {
995     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
996     return;
997   }
998
999   TString fileName("multiplicity");
1000   if (fSelectProcessType == 1)
1001     fileName += "ND";
1002   if (fSelectProcessType == 2)
1003     fileName += "SD";
1004   if (fSelectProcessType == 3)
1005     fileName += "DD";
1006   fileName += ".root";
1007   TFile* file = TFile::Open(fileName, "RECREATE");
1008
1009   fMultiplicity->SaveHistograms();
1010   for (Int_t i = 0; i < 8; ++i)
1011     if (fParticleCorrection[i])
1012       fParticleCorrection[i]->SaveHistograms();
1013   if (fParticleSpecies)
1014     fParticleSpecies->Write();
1015   if (fdNdpT)
1016     fdNdpT->Write();
1017
1018   fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1019   if (fTemp1)
1020     fTemp1->Write();
1021     
1022   fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
1023   if (fEtaPhi)
1024     fEtaPhi->Write();
1025   
1026   fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
1027   if (fVertex)
1028     fVertex->Write();
1029   
1030   for (Int_t i=0; i<3; i++)
1031   {
1032     fEta[i] = dynamic_cast<TH1*> (fOutput->FindObject(Form("fEta_%d", i)));
1033     if (fEta[i])
1034     {
1035       fEta[i]->Sumw2();
1036       Float_t events = fMultiplicity->GetMultiplicityESD(i)->Integral(1, fMultiplicity->GetMultiplicityESD(i)->GetNbinsX());
1037       if (events > 0)
1038         fEta[i]->Scale(1.0 / events);
1039       fEta[i]->Scale(1.0 / fEta[i]->GetXaxis()->GetBinWidth(1));
1040       fEta[i]->Write();
1041     }
1042   }
1043   
1044   TObjString option(fOption);
1045   option.Write();
1046
1047   file->Close();
1048
1049   Printf("Written result to %s", fileName.Data());
1050 }