adding flag for non-field data
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaTask.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaTask.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 <TH1F.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 #include <TGraph.h>
20
21 #include <AliLog.h>
22 #include <AliESDVertex.h>
23 #include <AliESDEvent.h>
24 #include <AliStack.h>
25 #include <AliHeader.h>
26 #include <AliGenEventHeader.h>
27 #include <AliMultiplicity.h>
28 #include <AliAnalysisManager.h>
29 #include <AliMCEventHandler.h>
30 #include <AliMCEvent.h>
31 #include <AliESDInputHandler.h>
32 #include <AliESDHeader.h>
33
34 #include "AliESDtrackCuts.h"
35 #include "AliPWG0Helper.h"
36 #include "AliCorrection.h"
37 #include "AliCorrectionMatrix3D.h"
38 #include "dNdEta/dNdEtaAnalysis.h"
39
40 ClassImp(AlidNdEtaTask)
41
42 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
43   AliAnalysisTask("AlidNdEtaTask", ""),
44   fESD(0),
45   fOutput(0),
46   fOption(opt),
47   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
48   fTrigger(AliPWG0Helper::kMB1),
49   fFillPhi(kFALSE),
50   fDeltaPhiCut(-1),
51   fReadMC(kFALSE),
52   fUseMCVertex(kFALSE),
53   fOnlyPrimaries(kFALSE),
54   fUseMCKine(kFALSE),
55   fEsdTrackCuts(0),
56   fdNdEtaAnalysisESD(0),
57   fMult(0),
58   fMultVtx(0),
59   fEvents(0),
60   fVertexResolution(0),
61   fdNdEtaAnalysis(0),
62   fdNdEtaAnalysisND(0),
63   fdNdEtaAnalysisNSD(0),
64   fdNdEtaAnalysisTr(0),
65   fdNdEtaAnalysisTrVtx(0),
66   fdNdEtaAnalysisTracks(0),
67   fPartPt(0),
68   fVertex(0),
69   fPhi(0),
70   fRawPt(0),
71   fEtaPhi(0),
72   fDeltaPhi(0),
73   fDeltaTheta(0),
74   fFiredChips(0),
75   fTriggerVsTime(0),
76   fStats(0)
77 {
78   //
79   // Constructor. Initialization of pointers
80   //
81
82   // Define input and output slots here
83   DefineInput(0, TChain::Class());
84   DefineOutput(0, TList::Class());
85   
86   fZPhi[0] = 0;
87   fZPhi[1] = 0;
88
89   AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
90 }
91
92 AlidNdEtaTask::~AlidNdEtaTask()
93 {
94   //
95   // Destructor
96   //
97
98   // histograms are in the output list and deleted when the output
99   // list is deleted by the TSelector dtor
100
101   if (fOutput) {
102     delete fOutput;
103     fOutput = 0;
104   }
105 }
106
107 Bool_t AlidNdEtaTask::Notify()
108 {
109   static Int_t count = 0;
110   count++;
111   Printf("Processing %d. file: %s", count, ((TTree*) GetInputData(0))->GetCurrentFile()->GetName());
112   return kTRUE;
113 }
114
115 //________________________________________________________________________
116 void AlidNdEtaTask::ConnectInputData(Option_t *)
117 {
118   // Connect ESD
119   // Called once
120
121   Printf("AlidNdEtaTask::ConnectInputData called");
122
123   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
124
125   if (!esdH) {
126     Printf("ERROR: Could not get ESDInputHandler");
127   } else {
128     fESD = esdH->GetEvent();
129     
130     TString branches("AliESDHeader Vertex ");
131
132     if (fAnalysisMode & AliPWG0Helper::kSPD || fTrigger == AliPWG0Helper::kOfflineMB1 || fTrigger == AliPWG0Helper::kOfflineMB2 || fTrigger == AliPWG0Helper::kOfflineMB3 || fTrigger == AliPWG0Helper::kOfflineFASTOR)
133       branches += "AliMultiplicity ";
134       
135     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
136       branches += "Tracks ";
137   
138     // Enable only the needed branches
139     esdH->SetActiveBranches(branches);
140   }
141
142   // disable info messages of AliMCEvent (per event)
143   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
144 }
145
146 void AlidNdEtaTask::CreateOutputObjects()
147 {
148   // create result objects and add to output list
149
150   Printf("AlidNdEtaTask::CreateOutputObjects");
151
152   if (fOnlyPrimaries)
153     Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
154
155   if (fUseMCKine)
156     Printf("WARNING: Using MC kine information. This is for systematical checks only.");
157
158   if (fUseMCVertex)
159     Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
160
161   fOutput = new TList;
162   fOutput->SetOwner();
163
164   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
165   fOutput->Add(fdNdEtaAnalysisESD);
166
167   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
168   fOutput->Add(fMult);
169
170   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
171   fOutput->Add(fMultVtx);
172
173   for (Int_t i=0; i<3; ++i)
174   {
175     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
176     fPartEta[i]->Sumw2();
177     fOutput->Add(fPartEta[i]);
178   }
179
180   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
181   fOutput->Add(fEvents);
182
183   Float_t resMax = (fAnalysisMode & AliPWG0Helper::kSPD) ? 0.2 : 2;
184   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, resMax);
185   fOutput->Add(fVertexResolution);
186
187   fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
188   fOutput->Add(fPhi);
189
190   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*5, 0, 2 * TMath::Pi());
191   fOutput->Add(fEtaPhi);
192
193   fTriggerVsTime = new TGraph; //TH1F("fTriggerVsTime", "fTriggerVsTime;trigger time;count", 100, 0, 100);
194   fTriggerVsTime->SetName("fTriggerVsTime");
195   fTriggerVsTime->GetXaxis()->SetTitle("trigger time");
196   fTriggerVsTime->GetYaxis()->SetTitle("count");
197   fOutput->Add(fTriggerVsTime);
198
199   fStats = new TH1F("fStats", "fStats", 3, 0.5, 3.5);
200   fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
201   fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
202   fStats->GetXaxis()->SetBinLabel(3, "trigger");
203   fOutput->Add(fStats);
204
205   if (fAnalysisMode & AliPWG0Helper::kSPD)
206   {
207     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.2, 0.2);
208     fOutput->Add(fDeltaPhi);
209     fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.2, 0.2);
210     fOutput->Add(fDeltaTheta);
211     fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
212     fOutput->Add(fFiredChips);
213     for (Int_t i=0; i<2; i++)
214     {
215       fZPhi[i] = new TH2F(Form("fZPhi_%d", i), Form("fZPhi Layer %d;z (cm);#phi (rad.)", i), 200, -20, 20, 180, 0, TMath::Pi() * 2);
216       fOutput->Add(fZPhi[i]);
217     }
218   }
219
220   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
221   {
222     fRawPt =  new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
223     fOutput->Add(fRawPt);
224   }
225
226   fVertex = new TH3F("vertex_check", "vertex_check", 100, -1, 1, 100, -1, 1, 100, -30, 30);
227   fOutput->Add(fVertex);
228
229   if (fReadMC)
230   {
231     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
232     fOutput->Add(fdNdEtaAnalysis);
233
234     fdNdEtaAnalysisND = new dNdEtaAnalysis("dndetaND", "dndetaND", fAnalysisMode);
235     fOutput->Add(fdNdEtaAnalysisND);
236
237     fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
238     fOutput->Add(fdNdEtaAnalysisNSD);
239
240     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
241     fOutput->Add(fdNdEtaAnalysisTr);
242
243     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
244     fOutput->Add(fdNdEtaAnalysisTrVtx);
245
246     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
247     fOutput->Add(fdNdEtaAnalysisTracks);
248
249     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
250     fPartPt->Sumw2();
251     fOutput->Add(fPartPt);
252   }
253
254   if (fEsdTrackCuts)
255   {
256     fEsdTrackCuts->SetName("fEsdTrackCuts");
257     fOutput->Add(fEsdTrackCuts);
258   }
259 }
260
261 void AlidNdEtaTask::Exec(Option_t*)
262 {
263   // process the event
264
265   // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
266   Bool_t eventTriggered = kFALSE;
267   const AliESDVertex* vtxESD = 0;
268
269   // post the data already here
270   PostData(0, fOutput);
271
272   // ESD analysis
273   if (fESD)
274   {
275     // check event type (should be PHYSICS = 7)
276     AliESDHeader* esdHeader = fESD->GetHeader();
277     if (!esdHeader)
278     {
279       Printf("ERROR: esdHeader could not be retrieved");
280       return;
281     }
282
283     /*
284     UInt_t eventType = esdHeader->GetEventType();
285     if (eventType != 7)
286     {
287       Printf("Skipping event because it is of type %d", eventType);
288       return;
289     }
290     */
291
292     // trigger definition
293     eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
294     if (eventTriggered)
295       fStats->Fill(3);
296
297     // get the ESD vertex
298     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
299
300     Double_t vtx[3];
301
302     // fill z vertex resolution
303     if (vtxESD)
304     {
305       fVertexResolution->Fill(vtxESD->GetZRes());
306       fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
307
308       if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
309       {
310           vtxESD->GetXYZ(vtx);
311
312           // vertex stats
313           if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
314           {
315             fStats->Fill(1);
316           }
317           else if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
318             fStats->Fill(2);
319       }
320       else
321         vtxESD = 0;
322     }
323
324     // needed for syst. studies
325     AliStack* stack = 0;
326     TArrayF vtxMC(3);
327
328     if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
329       if (!fReadMC) {
330         Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
331         return;
332       }
333
334       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
335       if (!eventHandler) {
336         Printf("ERROR: Could not retrieve MC event handler");
337         return;
338       }
339
340       AliMCEvent* mcEvent = eventHandler->MCEvent();
341       if (!mcEvent) {
342         Printf("ERROR: Could not retrieve MC event");
343         return;
344       }
345
346       AliHeader* header = mcEvent->Header();
347       if (!header)
348       {
349         AliDebug(AliLog::kError, "Header not available");
350         return;
351       }
352
353       // get the MC vertex
354       AliGenEventHeader* genHeader = header->GenEventHeader();
355       if (!genHeader)
356       {
357         AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
358         return;
359       }
360       genHeader->PrimaryVertex(vtxMC);
361
362       if (fUseMCVertex)
363         vtx[2] = vtxMC[2];
364
365       stack = mcEvent->Stack();
366       if (!stack)
367       {
368         AliDebug(AliLog::kError, "Stack not available");
369         return;
370       }
371     }
372
373     // create list of (label, eta, pt) tuples
374     Int_t inputCount = 0;
375     Int_t* labelArr = 0;
376     Float_t* etaArr = 0;
377     Float_t* thirdDimArr = 0;
378     if (fAnalysisMode & AliPWG0Helper::kSPD)
379     {
380       // get tracklets
381       const AliMultiplicity* mult = fESD->GetMultiplicity();
382       if (!mult)
383       {
384         AliDebug(AliLog::kError, "AliMultiplicity not available");
385         return;
386       }
387
388       labelArr = new Int_t[mult->GetNumberOfTracklets()];
389       etaArr = new Float_t[mult->GetNumberOfTracklets()];
390       thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
391
392       // get multiplicity from SPD tracklets
393       for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
394       {
395         //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
396
397         if (fOnlyPrimaries)
398           if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
399             continue;
400
401         Float_t deltaPhi = mult->GetDeltaPhi(i);
402         // prevent values to be shifted by 2 Pi()
403         if (deltaPhi < -TMath::Pi())
404           deltaPhi += TMath::Pi() * 2;
405         if (deltaPhi > TMath::Pi())
406           deltaPhi -= TMath::Pi() * 2;
407
408         if (TMath::Abs(deltaPhi) > 1)
409           printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
410
411         Int_t label = mult->GetLabel(i, 0);
412         Float_t eta = mult->GetEta(i);
413
414         // control histograms
415         Float_t phi = mult->GetPhi(i);
416         if (phi < 0)
417           phi += TMath::Pi() * 2;
418         fPhi->Fill(phi);
419         fEtaPhi->Fill(eta, phi);
420         
421         if (deltaPhi < 0.01)
422         {
423           // layer 0
424           Float_t z = vtx[2] + 3.9 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta)));
425           fZPhi[0]->Fill(z, phi);
426           // layer 1
427           z = vtx[2] + 7.6 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta)));
428           fZPhi[1]->Fill(z, phi);
429         }
430
431         fDeltaPhi->Fill(deltaPhi);
432         fDeltaTheta->Fill(mult->GetDeltaTheta(i));
433
434         if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
435           continue;
436
437         if (fUseMCKine)
438         {
439           if (label > 0)
440           {
441             TParticle* particle = stack->Particle(label);
442             eta = particle->Eta();
443             phi = particle->Phi();
444           }
445           else
446             Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
447         }
448         
449         etaArr[inputCount] = eta;
450         labelArr[inputCount] = label;
451         thirdDimArr[inputCount] = phi;
452         ++inputCount;
453       }
454
455       if (!fFillPhi)
456       {
457         // fill multiplicity in third axis
458         for (Int_t i=0; i<inputCount; ++i)
459           thirdDimArr[i] = inputCount;
460       }
461
462       Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
463       fFiredChips->Fill(firedChips, inputCount);
464       Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
465     }
466     else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
467     {
468       if (!fEsdTrackCuts)
469       {
470         AliDebug(AliLog::kError, "fESDTrackCuts not available");
471         return;
472       }
473
474       if (vtxESD)
475       {
476         // get multiplicity from ESD tracks
477         TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
478         Int_t nGoodTracks = list->GetEntries();
479         Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
480   
481         labelArr = new Int_t[nGoodTracks];
482         etaArr = new Float_t[nGoodTracks];
483         thirdDimArr = new Float_t[nGoodTracks];
484   
485         // loop over esd tracks
486         for (Int_t i=0; i<nGoodTracks; i++)
487         {
488           AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
489           if (!esdTrack)
490           {
491             AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
492             continue;
493           }
494           
495           Float_t phi = esdTrack->Phi();
496           if (phi < 0)
497             phi += TMath::Pi() * 2;
498   
499           Float_t eta = esdTrack->Eta();
500           Int_t label = TMath::Abs(esdTrack->GetLabel());
501           Float_t pT  = esdTrack->Pt();
502           
503           // force pT to fixed value without B field
504           if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
505             pT = 1;
506   
507           fPhi->Fill(phi);
508           fEtaPhi->Fill(eta, phi);
509           if (eventTriggered && vtxESD)
510             fRawPt->Fill(pT);
511   
512           if (fOnlyPrimaries)
513           {
514             if (label == 0)
515               continue;
516             
517             if (stack->IsPhysicalPrimary(label) == kFALSE)
518               continue;
519           }
520   
521           if (fUseMCKine)
522           {
523             if (label > 0)
524             {
525               TParticle* particle = stack->Particle(label);
526               eta = particle->Eta();
527               pT = particle->Pt();
528             }
529             else
530               Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
531           }
532   
533           etaArr[inputCount] = eta;
534           labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
535           thirdDimArr[inputCount] = pT;
536           ++inputCount;
537         }
538         
539         // TODO restrict inputCount used as measure for the multiplicity to |eta| < 1
540   
541         delete list;
542       }
543     }
544     else
545       return;
546
547     // Processing of ESD information (always)
548     if (eventTriggered)
549     {
550       // control hist
551       fMult->Fill(inputCount);
552       fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
553
554       fTriggerVsTime->SetPoint(fTriggerVsTime->GetN(), fESD->GetTimeStamp(), 1);
555
556       if (vtxESD)
557       {
558         // control hist
559         fMultVtx->Fill(inputCount);
560
561         for (Int_t i=0; i<inputCount; ++i)
562         {
563           Float_t eta = etaArr[i];
564           Float_t thirdDim = thirdDimArr[i];
565
566           fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, thirdDim);
567
568           if (TMath::Abs(vtx[2]) < 10)
569           {
570             fPartEta[0]->Fill(eta);
571
572             if (vtx[2] < 0)
573               fPartEta[1]->Fill(eta);
574             else
575               fPartEta[2]->Fill(eta);
576           }
577         }
578
579         // for event count per vertex
580         fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
581
582         // control hist
583         fEvents->Fill(vtx[2]);
584
585         if (fReadMC)
586         {
587           // from tracks is only done for triggered and vertex reconstructed events
588           for (Int_t i=0; i<inputCount; ++i)
589           {
590             Int_t label = labelArr[i];
591
592             if (label < 0)
593               continue;
594
595             //Printf("Getting particle of track %d", label);
596             TParticle* particle = stack->Particle(label);
597
598             if (!particle)
599             {
600               AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
601               continue;
602             }
603
604             Float_t thirdDim = -1;
605             if (fAnalysisMode & AliPWG0Helper::kSPD)
606             {
607               if (fFillPhi)
608               {
609                 thirdDim = particle->Phi();
610               }
611               else
612                 thirdDim = inputCount;
613             }
614             else
615               thirdDim = particle->Pt();
616
617             fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
618           } // end of track loop
619
620           // for event count per vertex
621           fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
622         }
623       }
624     }
625
626     if (etaArr)
627       delete[] etaArr;
628     if (labelArr)
629       delete[] labelArr;
630     if (thirdDimArr)
631       delete[] thirdDimArr;
632   }
633
634   if (fReadMC)   // Processing of MC information (optional)
635   {
636     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
637     if (!eventHandler) {
638       Printf("ERROR: Could not retrieve MC event handler");
639       return;
640     }
641
642     AliMCEvent* mcEvent = eventHandler->MCEvent();
643     if (!mcEvent) {
644       Printf("ERROR: Could not retrieve MC event");
645       return;
646     }
647
648     AliStack* stack = mcEvent->Stack();
649     if (!stack)
650     {
651       AliDebug(AliLog::kError, "Stack not available");
652       return;
653     }
654
655     AliHeader* header = mcEvent->Header();
656     if (!header)
657     {
658       AliDebug(AliLog::kError, "Header not available");
659       return;
660     }
661
662     // get the MC vertex
663     AliGenEventHeader* genHeader = header->GenEventHeader();
664     if (!genHeader)
665     {
666       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
667       return;
668     }
669
670     TArrayF vtxMC(3);
671     genHeader->PrimaryVertex(vtxMC);
672
673     // get process type
674     Int_t processType = AliPWG0Helper::GetEventProcessType(header);
675     AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
676
677     if (processType==AliPWG0Helper::kInvalidProcess)
678       AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
679
680     // loop over mc particles
681     Int_t nPrim  = stack->GetNprimary();
682
683     Int_t nAcceptedParticles = 0;
684
685     // count particles first, then fill
686     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
687     {
688       if (!stack->IsPhysicalPrimary(iMc))
689         continue;
690
691       //Printf("Getting particle %d", iMc);
692       TParticle* particle = stack->Particle(iMc);
693
694       if (!particle)
695         continue;
696
697       if (particle->GetPDG()->Charge() == 0)
698         continue;
699
700       //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
701       Float_t eta = particle->Eta();
702
703        // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
704       if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
705         nAcceptedParticles++;
706     }
707
708     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
709     {
710       if (!stack->IsPhysicalPrimary(iMc))
711         continue;
712
713       //Printf("Getting particle %d", iMc);
714       TParticle* particle = stack->Particle(iMc);
715
716       if (!particle)
717         continue;
718
719       if (particle->GetPDG()->Charge() == 0)
720         continue;
721
722       Float_t eta = particle->Eta();
723       Float_t thirdDim = -1;
724
725       if (fAnalysisMode & AliPWG0Helper::kSPD)
726       {
727         if (fFillPhi)
728         {
729           thirdDim = particle->Phi();
730         }
731         else
732           thirdDim = nAcceptedParticles;
733       }
734       else
735         thirdDim = particle->Pt();
736
737       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
738
739       if (processType != AliPWG0Helper::kSD)
740         fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
741
742       if (processType == AliPWG0Helper::kND)
743         fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
744
745       if (eventTriggered)
746       {
747         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
748         if (vtxESD)
749           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
750       }
751
752       if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0)
753         fPartPt->Fill(particle->Pt());
754     }
755
756     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
757     if (processType != AliPWG0Helper::kSD)
758       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
759     if (processType == AliPWG0Helper::kND)
760       fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
761
762     if (eventTriggered)
763     {
764       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
765       if (vtxESD)
766         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
767     }
768   }
769 }
770
771 void AlidNdEtaTask::Terminate(Option_t *)
772 {
773   // The Terminate() function is the last function to be called during
774   // a query. It always runs on the client, it can be used to present
775   // the results graphically or save the results to file.
776
777   fOutput = dynamic_cast<TList*> (GetOutputData(0));
778   if (!fOutput)
779     Printf("ERROR: fOutput not available");
780
781   if (fOutput)
782   {
783     fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
784     fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
785     fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
786     for (Int_t i=0; i<3; ++i)
787       fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
788     fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
789     fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
790
791     fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
792     fRawPt = dynamic_cast<TH1F*> (fOutput->FindObject("fRawPt"));
793     fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
794     for (Int_t i=0; i<2; ++i)
795       fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
796     fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
797     fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
798     fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
799     fTriggerVsTime = dynamic_cast<TGraph*> (fOutput->FindObject("fTriggerVsTime"));
800     fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
801
802     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
803   }
804
805   if (!fdNdEtaAnalysisESD)
806   {
807     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
808   }
809   else
810   {
811     if (fMult && fMultVtx)
812     {
813       new TCanvas;
814       fMult->Draw();
815       fMultVtx->SetLineColor(2);
816       fMultVtx->Draw("SAME");
817     }
818
819     if (fFiredChips)
820     {
821       new TCanvas;
822       fFiredChips->Draw("COLZ");
823     }
824
825     if (fPartEta[0])
826     {
827       Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
828       Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
829
830       Printf("%d events with vertex used", events1 + events2);
831
832       if (events1 > 0 && events2 > 0)
833       {
834         fPartEta[0]->Scale(1.0 / (events1 + events2));
835         fPartEta[1]->Scale(1.0 / events1);
836         fPartEta[2]->Scale(1.0 / events2);
837
838         for (Int_t i=0; i<3; ++i)
839           fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
840
841         new TCanvas("control", "control", 500, 500);
842         for (Int_t i=0; i<3; ++i)
843         {
844           fPartEta[i]->SetLineColor(i+1);
845           fPartEta[i]->Draw((i==0) ? "" : "SAME");
846         }
847       }
848     }
849
850     if (fEvents)
851     {
852         new TCanvas("control3", "control3", 500, 500);
853         fEvents->Draw();
854     }
855
856     TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
857
858     if (fdNdEtaAnalysisESD)
859       fdNdEtaAnalysisESD->SaveHistograms();
860
861     if (fEsdTrackCuts)
862       fEsdTrackCuts->SaveHistograms("esd_track_cuts");
863
864     if (fMult)
865       fMult->Write();
866
867     if (fMultVtx)
868       fMultVtx->Write();
869
870     for (Int_t i=0; i<3; ++i)
871       if (fPartEta[i])
872         fPartEta[i]->Write();
873
874     if (fEvents)
875       fEvents->Write();
876
877     if (fVertexResolution)
878       fVertexResolution->Write();
879
880     if (fDeltaPhi)
881       fDeltaPhi->Write();
882
883     if (fDeltaTheta)
884       fDeltaTheta->Write();
885     
886     if (fPhi)
887       fPhi->Write();
888
889     if (fRawPt)
890       fRawPt->Write();
891
892     if (fEtaPhi)
893       fEtaPhi->Write();
894
895     for (Int_t i=0; i<2; ++i)
896       if (fZPhi[i])
897         fZPhi[i]->Write();
898     
899     if (fFiredChips)
900       fFiredChips->Write();
901
902     if (fTriggerVsTime)
903       fTriggerVsTime->Write();
904
905     if (fStats)
906       fStats->Write();
907
908     fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
909     if (fVertex)
910       fVertex->Write();
911
912     fout->Write();
913     fout->Close();
914
915     Printf("Writting result to analysis_esd_raw.root");
916   }
917
918   if (fReadMC)
919   {
920     if (fOutput)
921     {
922       fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
923       fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
924       fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
925       fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
926       fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
927       fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
928       fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
929     }
930
931     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
932     {
933       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
934       return;
935     }
936
937     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
938     fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
939     fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
940     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
941     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
942     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
943
944     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
945     fPartPt->Scale(1.0/events);
946     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
947
948     TFile* fout = new TFile("analysis_mc.root","RECREATE");
949
950     fdNdEtaAnalysis->SaveHistograms();
951     fdNdEtaAnalysisND->SaveHistograms();
952     fdNdEtaAnalysisNSD->SaveHistograms();
953     fdNdEtaAnalysisTr->SaveHistograms();
954     fdNdEtaAnalysisTrVtx->SaveHistograms();
955     fdNdEtaAnalysisTracks->SaveHistograms();
956
957     if (fPartPt)
958       fPartPt->Write();
959
960     fout->Write();
961     fout->Close();
962
963     Printf("Writting result to analysis_mc.root");
964
965     if (fPartPt)
966     {
967       new TCanvas("control2", "control2", 500, 500);
968       fPartPt->Draw();
969     }
970   }
971 }