updated dndeta analysis
[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 <AliESDInputHandlerRP.h>
33 #include <AliESDHeader.h>
34
35 #include "AliESDtrackCuts.h"
36 #include "AliPWG0Helper.h"
37 #include "AliCorrection.h"
38 #include "AliCorrectionMatrix3D.h"
39 #include "dNdEta/dNdEtaAnalysis.h"
40 #include "AliTriggerAnalysis.h"
41 #include "AliPhysicsSelection.h"
42
43 //#define FULLALIROOT
44
45 #ifdef FULLALIROOT
46   #include "../ITS/AliITSRecPoint.h"
47   #include "AliCDBManager.h"
48   #include "AliCDBEntry.h"
49   #include "AliGeomManager.h"
50   #include "TGeoManager.h"
51 #endif
52
53 ClassImp(AlidNdEtaTask)
54
55 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
56   AliAnalysisTask("AlidNdEtaTask", ""),
57   fESD(0),
58   fOutput(0),
59   fOption(opt),
60   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
61   fTrigger(AliTriggerAnalysis::kMB1),
62   fRequireTriggerClass(),
63   fRejectTriggerClass(),
64   fFillPhi(kFALSE),
65   fDeltaPhiCut(-1),
66   fReadMC(kFALSE),
67   fUseMCVertex(kFALSE),
68   fOnlyPrimaries(kFALSE),
69   fUseMCKine(kFALSE),
70   fCheckEventType(kFALSE),
71   fSymmetrize(kFALSE),
72   fMultAxisEta1(kFALSE),
73   fDiffTreatment(AliPWG0Helper::kMCFlags),
74   fEsdTrackCuts(0),
75   fTriggerAnalysis(0),
76   fdNdEtaAnalysisESD(0),
77   fMult(0),
78   fMultVtx(0),
79   fEvents(0),
80   fVertexResolution(0),
81   fdNdEtaAnalysis(0),
82   fdNdEtaAnalysisND(0),
83   fdNdEtaAnalysisNSD(0),
84   fdNdEtaAnalysisOnePart(0),
85   fdNdEtaAnalysisTr(0),
86   fdNdEtaAnalysisTrVtx(0),
87   fdNdEtaAnalysisTracks(0),
88   fPartPt(0),
89   fVertex(0),
90   fVertexVsMult(0),
91   fPhi(0),
92   fRawPt(0),
93   fEtaPhi(0),
94   fModuleMap(0),
95   fDeltaPhi(0),
96   fDeltaTheta(0),
97   fFiredChips(0),
98   fTrackletsVsClusters(0),
99   fTrackletsVsUnassigned(0),
100   fStats(0),
101   fStats2(0)
102 {
103   //
104   // Constructor. Initialization of pointers
105   //
106
107   // Define input and output slots here
108   DefineInput(0, TChain::Class());
109   DefineOutput(0, TList::Class());
110   
111   fZPhi[0] = 0;
112   fZPhi[1] = 0;
113
114   AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
115 }
116
117 AlidNdEtaTask::~AlidNdEtaTask()
118 {
119   //
120   // Destructor
121   //
122
123   // histograms are in the output list and deleted when the output
124   // list is deleted by the TSelector dtor
125
126   if (fOutput) {
127     delete fOutput;
128     fOutput = 0;
129   }
130 }
131
132 Bool_t AlidNdEtaTask::Notify()
133 {
134   static Int_t count = 0;
135   count++;
136   Printf("Processing %d. file: %s", count, ((TTree*) GetInputData(0))->GetCurrentFile()->GetName());
137   return kTRUE;
138 }
139
140 //________________________________________________________________________
141 void AlidNdEtaTask::ConnectInputData(Option_t *)
142 {
143   // Connect ESD
144   // Called once
145
146   Printf("AlidNdEtaTask::ConnectInputData called");
147
148   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
149
150   if (!esdH) {
151     Printf("ERROR: Could not get ESDInputHandler");
152   } else {
153     fESD = esdH->GetEvent();
154     
155     TString branches("AliESDHeader Vertex ");
156
157     if (fAnalysisMode & AliPWG0Helper::kSPD || fTrigger & AliTriggerAnalysis::kOfflineFlag)
158       branches += "AliMultiplicity ";
159       
160     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
161       branches += "Tracks ";
162   
163     // Enable only the needed branches
164     esdH->SetActiveBranches(branches);
165   }
166
167   // disable info messages of AliMCEvent (per event)
168   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
169   
170   #ifdef FULLALIROOT
171     if (fCheckEventType)
172       AliCDBManager::Instance()->SetDefaultStorage("raw://");
173     else
174       AliCDBManager::Instance()->SetDefaultStorage("MC", "Residual");
175     AliCDBManager::Instance()->SetRun(0);
176     
177     AliCDBManager* mgr = AliCDBManager::Instance();
178     AliCDBEntry* obj = mgr->Get(AliCDBPath("GRP", "Geometry", "Data"));
179     AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
180     
181     AliGeomManager::GetNalignable("ITS");
182     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
183   #endif
184 }
185
186 void AlidNdEtaTask::CreateOutputObjects()
187 {
188   // create result objects and add to output list
189
190   Printf("AlidNdEtaTask::CreateOutputObjects");
191
192   if (fOnlyPrimaries)
193     Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
194
195   if (fUseMCKine)
196     Printf("WARNING: Using MC kine information. This is for systematical checks only.");
197
198   if (fUseMCVertex)
199     Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
200
201   fOutput = new TList;
202   fOutput->SetOwner();
203
204   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
205   fOutput->Add(fdNdEtaAnalysisESD);
206
207   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
208   fOutput->Add(fMult);
209
210   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
211   fOutput->Add(fMultVtx);
212
213   for (Int_t i=0; i<3; ++i)
214   {
215     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -3, 3);
216     fPartEta[i]->Sumw2();
217     fOutput->Add(fPartEta[i]);
218   }
219
220   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
221   fOutput->Add(fEvents);
222
223   fVertexResolution = new TH2F("dndeta_vertex_resolution_z", ";z resolution;#Delta #phi", 1000, 0, 2, 100, 0, 0.2);
224   fOutput->Add(fVertexResolution);
225
226   fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
227   fOutput->Add(fPhi);
228
229   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
230   fOutput->Add(fEtaPhi);
231   
232   fStats = new TH1F("fStats", "fStats", 5, 0.5, 5.5);
233   fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
234   fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
235   fStats->GetXaxis()->SetBinLabel(3, "trigger");
236   fStats->GetXaxis()->SetBinLabel(4, "physics events");
237   fStats->GetXaxis()->SetBinLabel(5, "physics events after veto");
238   fOutput->Add(fStats);
239   
240   fStats2 = new TH2F("fStats2", "fStats2", 7, -0.5, 6.5, 10, -0.5, 9.5);
241   
242   fStats2->GetXaxis()->SetBinLabel(1, "No trigger");
243   fStats2->GetXaxis()->SetBinLabel(2, "No Vertex");
244   fStats2->GetXaxis()->SetBinLabel(3, "Vtx res. too bad");
245   fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 15");
246   fStats2->GetXaxis()->SetBinLabel(5, "|z-vtx| > 10");
247   fStats2->GetXaxis()->SetBinLabel(6, "0 tracklets");
248   fStats2->GetXaxis()->SetBinLabel(7, "Selected");
249   
250   fStats2->GetYaxis()->SetBinLabel(1, "n/a");
251   fStats2->GetYaxis()->SetBinLabel(2, "empty");
252   fStats2->GetYaxis()->SetBinLabel(3, "BBA");
253   fStats2->GetYaxis()->SetBinLabel(4, "BBC");
254   fStats2->GetYaxis()->SetBinLabel(5, "BBA BBC");
255   fStats2->GetYaxis()->SetBinLabel(6, "BGA");
256   fStats2->GetYaxis()->SetBinLabel(7, "BGC");
257   fStats2->GetYaxis()->SetBinLabel(8, "BGA BGC");
258   fStats2->GetYaxis()->SetBinLabel(9, "BBA BGC");
259   fStats2->GetYaxis()->SetBinLabel(10, "BGA BBC");
260   fOutput->Add(fStats2);
261
262   fTrackletsVsClusters = new TH2F("fTrackletsVsClusters", ";tracklets;clusters in ITS", 50, -0.5, 49.5, 1000, -0.5, 999.5);
263   fOutput->Add(fTrackletsVsClusters);
264   
265   if (fAnalysisMode & AliPWG0Helper::kSPD)
266   {
267     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.16, 0.16);
268     fOutput->Add(fDeltaPhi);
269     fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.05, 0.05);
270     fOutput->Add(fDeltaTheta);
271     fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
272     fOutput->Add(fFiredChips);
273     fTrackletsVsUnassigned = new TH2F("fTrackletsVsUnassigned", ";tracklets;unassigned clusters in L0", 50, -0.5, 49.5, 200, -0.5, 199.5);
274     fOutput->Add(fTrackletsVsUnassigned);
275     for (Int_t i=0; i<2; i++)
276     {
277       fZPhi[i] = new TH2F(Form("fZPhi_%d", i), Form("fZPhi Layer %d;z (cm);#phi (rad.)", i), 200, -15, 15, 200, 0, TMath::Pi() * 2);
278       fOutput->Add(fZPhi[i]);
279     }
280     
281     fModuleMap = new TH1F("fModuleMap", "fModuleMap;module number;cluster count", 240, -0.5, 239.5);
282     fOutput->Add(fModuleMap);
283   }
284
285   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
286   {
287     fRawPt =  new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
288     fOutput->Add(fRawPt);
289   }
290
291   fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
292   fOutput->Add(fVertex);
293   
294   fVertexVsMult = new TH3F("fVertexVsMult", "fVertexVsMult;x;y;multiplicity", 100, -1, 1, 100, -1, 1, 100, -0.5, 99.5);
295   fOutput->Add(fVertexVsMult);
296
297   if (fReadMC)
298   {
299     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
300     fOutput->Add(fdNdEtaAnalysis);
301
302     fdNdEtaAnalysisND = new dNdEtaAnalysis("dndetaND", "dndetaND", fAnalysisMode);
303     fOutput->Add(fdNdEtaAnalysisND);
304
305     fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
306     fOutput->Add(fdNdEtaAnalysisNSD);
307
308     fdNdEtaAnalysisOnePart = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart", fAnalysisMode);
309     fOutput->Add(fdNdEtaAnalysisOnePart);
310
311     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
312     fOutput->Add(fdNdEtaAnalysisTr);
313
314     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
315     fOutput->Add(fdNdEtaAnalysisTrVtx);
316
317     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
318     fOutput->Add(fdNdEtaAnalysisTracks);
319
320     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
321     fPartPt->Sumw2();
322     fOutput->Add(fPartPt);
323   }
324
325   if (fEsdTrackCuts)
326   {
327     fEsdTrackCuts->SetName("fEsdTrackCuts");
328     fOutput->Add(fEsdTrackCuts);
329   }
330   
331   //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
332 }
333
334 void AlidNdEtaTask::Exec(Option_t*)
335 {
336   // process the event
337
338   // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
339   Bool_t eventTriggered = kFALSE;
340   const AliESDVertex* vtxESD = 0;
341
342   // post the data already here
343   PostData(0, fOutput);
344
345   // ESD analysis
346   if (fESD)
347   {
348     AliESDHeader* esdHeader = fESD->GetHeader();
349     if (!esdHeader)
350     {
351       Printf("ERROR: esdHeader could not be retrieved");
352       return;
353     }
354     
355     AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
356     if (!inputHandler)
357     {
358       Printf("ERROR: Could not receive input handler");
359       return;
360     }
361     
362     eventTriggered = inputHandler->IsEventSelected();
363         
364     if (!fTriggerAnalysis)
365     {
366       AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
367       if (physicsSelection)
368         fTriggerAnalysis = physicsSelection->GetTriggerAnalysis();
369     }
370       
371     if (eventTriggered)
372       eventTriggered = fTriggerAnalysis->IsTriggerFired(fESD, fTrigger);
373     
374     AliTriggerAnalysis::V0Decision v0A = fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kFALSE);
375     AliTriggerAnalysis::V0Decision v0C = fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kFALSE);
376     
377     Int_t vZero = 0;
378     if (v0A != AliTriggerAnalysis::kV0Invalid && v0C != AliTriggerAnalysis::kV0Invalid)
379     {
380       if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0Empty) vZero = 1;
381       if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0Empty) vZero = 2;
382       if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BB)    vZero = 3;
383       if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0BB)    vZero = 4;
384       if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0Empty) vZero = 5;
385       if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BG)    vZero = 6;
386       if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0BG)    vZero = 7;
387       if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0BG)    vZero = 8;
388       if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0BB)    vZero = 9;
389     }
390
391     if (vZero == 0)
392       Printf("VZERO: %d %d %d %d", vZero, eventTriggered, v0A, v0C);
393       
394     Bool_t filled = kFALSE;
395       
396     if (!eventTriggered)
397     {
398       fStats2->Fill(0.0, vZero);
399       filled = kTRUE;
400     }
401     
402     if (eventTriggered)
403       fStats->Fill(3);
404       
405     /*
406     // ITS cluster tree
407     AliESDInputHandlerRP* handlerRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalys
408 isManager()->GetInputEventHandler());
409     if (handlerRP)
410     {
411       TTree* itsClusterTree = handlerRP->GetTreeR("ITS");
412       if (!itsClusterTree)
413         return;
414
415       TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
416       TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
417
418       itsClusterBranch->SetAddress(&itsClusters);
419
420       Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
421
422       Int_t totalClusters = 0;
423
424       // loop over the its subdetectors
425       for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
426
427         if (!itsClusterTree->GetEvent(iIts))
428           continue;
429
430         Int_t nClusters = itsClusters->GetEntriesFast();
431         totalClusters += nClusters;
432
433         #ifdef FULLALIROOT
434           if (fAnalysisMode & AliPWG0Helper::kSPD)
435           {
436             // loop over clusters
437             while (nClusters--) {
438               AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
439
440               Int_t layer = cluster->GetLayer();
441
442               if (layer > 1)
443                 continue;
444
445               Float_t xyz[3] = {0., 0., 0.};
446               cluster->GetGlobalXYZ(xyz);
447
448               Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]);
449               Float_t z = xyz[2];
450
451               fZPhi[layer]->Fill(z, phi);
452               fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex());
453             }
454           }
455         #endif
456       }
457     }
458     */
459       
460     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
461     if (!vtxESD)
462     {
463       if (!filled)
464       {
465         fStats2->Fill(1, vZero);
466         filled = kTRUE;
467       }
468     }
469     else
470     {
471       if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
472       {
473         if (!filled)
474         {
475           fStats2->Fill(2, vZero);
476           filled = kTRUE;
477         }
478       }
479       
480       Double_t vtx[3];
481       vtxESD->GetXYZ(vtx);
482       
483       // try to compare spd vertex and vertexer from tracks
484       // remove vertices outside +- 15 cm
485       if (TMath::Abs(vtx[2]) > 15)
486       {
487         if (!filled)
488         {
489           fStats2->Fill(3, vZero);
490           filled = kTRUE;
491         }
492       }
493       
494       if (TMath::Abs(vtx[2]) > 10)
495       {
496         if (!filled)
497         {
498           fStats2->Fill(4, vZero);
499           filled = kTRUE;
500         }
501       }
502         
503       const AliMultiplicity* mult = fESD->GetMultiplicity();
504       if (!mult)
505         return;
506       
507       if (mult->GetNumberOfTracklets() == 0)
508       {
509         if (!filled)
510         {
511           fStats2->Fill(5, vZero);
512           filled = kTRUE;
513         }
514       }
515     }
516     
517     if (!filled)
518     {
519       fStats2->Fill(6, vZero);
520       //Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
521     }
522       
523     if (eventTriggered)
524       fStats->Fill(3);
525     
526     fStats->Fill(5);
527     
528     // get the ESD vertex
529     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
530     
531     Double_t vtx[3];
532
533     // fill z vertex resolution
534     if (vtxESD)
535     {
536       if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
537         fVertexResolution->Fill(vtxESD->GetZRes(), vtxESD->GetDispersion());
538       else
539         fVertexResolution->Fill(vtxESD->GetZRes(), 0);
540       
541       //if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
542       {
543         fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
544       }
545       
546       if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
547       {
548           vtxESD->GetXYZ(vtx);
549
550           // vertex stats
551           if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
552           {
553             fStats->Fill(1);
554             if (fCheckEventType && TMath::Abs(vtx[0]) > 0.3)
555             {
556               Printf("Suspicious x-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
557             }
558             if (fCheckEventType && (vtx[1] < 0.05 || vtx[1] > 0.5))
559             {
560               Printf("Suspicious y-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
561             }
562           }
563           
564           if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
565             fStats->Fill(2);
566           
567           // remove vertices outside +- 15 cm
568           if (TMath::Abs(vtx[2]) > 15)
569             vtxESD = 0;
570       }
571       else
572         vtxESD = 0;
573     }
574
575     // needed for syst. studies
576     AliStack* stack = 0;
577     TArrayF vtxMC(3);
578
579     if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
580       if (!fReadMC) {
581         Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
582         return;
583       }
584
585       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
586       if (!eventHandler) {
587         Printf("ERROR: Could not retrieve MC event handler");
588         return;
589       }
590
591       AliMCEvent* mcEvent = eventHandler->MCEvent();
592       if (!mcEvent) {
593         Printf("ERROR: Could not retrieve MC event");
594         return;
595       }
596
597       AliHeader* header = mcEvent->Header();
598       if (!header)
599       {
600         AliDebug(AliLog::kError, "Header not available");
601         return;
602       }
603
604       // get the MC vertex
605       AliGenEventHeader* genHeader = header->GenEventHeader();
606       if (!genHeader)
607       {
608         AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
609         return;
610       }
611       genHeader->PrimaryVertex(vtxMC);
612
613       if (fUseMCVertex)
614         vtx[2] = vtxMC[2];
615
616       stack = mcEvent->Stack();
617       if (!stack)
618       {
619         AliDebug(AliLog::kError, "Stack not available");
620         return;
621       }
622     }
623     
624     // create list of (label, eta, pt) tuples
625     Int_t inputCount = 0;
626     Int_t* labelArr = 0;
627     Float_t* etaArr = 0;
628     Float_t* thirdDimArr = 0;
629     if (fAnalysisMode & AliPWG0Helper::kSPD)
630     {
631       if (vtxESD)
632       {
633         // get tracklets
634         const AliMultiplicity* mult = fESD->GetMultiplicity();
635         if (!mult)
636         {
637           AliDebug(AliLog::kError, "AliMultiplicity not available");
638           return;
639         }
640   
641         Int_t arrayLength = mult->GetNumberOfTracklets();
642         if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
643           arrayLength += mult->GetNumberOfSingleClusters();
644           
645         labelArr = new Int_t[arrayLength];
646         etaArr = new Float_t[arrayLength];
647         thirdDimArr = new Float_t[arrayLength];
648   
649         // get multiplicity from SPD tracklets
650         for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
651         {
652           //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
653   
654           if (fOnlyPrimaries)
655             if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
656               continue;
657   
658           Float_t deltaPhi = mult->GetDeltaPhi(i);
659           // prevent values to be shifted by 2 Pi()
660           if (deltaPhi < -TMath::Pi())
661             deltaPhi += TMath::Pi() * 2;
662           if (deltaPhi > TMath::Pi())
663             deltaPhi -= TMath::Pi() * 2;
664   
665           if (TMath::Abs(deltaPhi) > 1)
666             printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
667   
668           Int_t label = mult->GetLabel(i, 0);
669           Float_t eta = mult->GetEta(i);
670           
671           // control histograms
672           Float_t phi = mult->GetPhi(i);
673           if (phi < 0)
674             phi += TMath::Pi() * 2;
675             
676           // TEST exclude probably inefficient phi region
677           //if (phi > 5.70 || phi < 0.06)
678           //  continue;
679             
680           fPhi->Fill(phi);
681           
682           if (vtxESD && TMath::Abs(vtx[2]) < 10)
683           {
684             fEtaPhi->Fill(eta, phi);
685             fDeltaPhi->Fill(deltaPhi);
686             fDeltaTheta->Fill(mult->GetDeltaTheta(i));
687           }
688           
689           if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
690             continue;
691   
692           if (fUseMCKine)
693           {
694             if (label > 0)
695             {
696               TParticle* particle = stack->Particle(label);
697               eta = particle->Eta();
698               phi = particle->Phi();
699             }
700             else
701               Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
702           }
703           
704           if (fSymmetrize)
705             eta = TMath::Abs(eta);
706   
707           etaArr[inputCount] = eta;
708           labelArr[inputCount] = label;
709           thirdDimArr[inputCount] = phi;
710           ++inputCount;
711         }
712         
713         if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
714         {
715           // get additional clusters from L0 
716           for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
717           {
718             etaArr[inputCount] = -TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.));
719             labelArr[inputCount] = -1;
720             thirdDimArr[inputCount] = mult->GetPhiSingle(i);
721             
722             ++inputCount;
723           }
724         }
725         
726         if (!fFillPhi)
727         {
728           // fill multiplicity in third axis
729           for (Int_t i=0; i<inputCount; ++i)
730             thirdDimArr[i] = inputCount;
731         }
732   
733         Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
734         fFiredChips->Fill(firedChips, inputCount);
735         Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
736         
737         fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters());
738       }
739     }
740     else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
741     {
742       if (!fEsdTrackCuts)
743       {
744         AliDebug(AliLog::kError, "fESDTrackCuts not available");
745         return;
746       }
747
748       Bool_t foundInEta10 = kFALSE;
749       
750       if (vtxESD)
751       {
752         // get multiplicity from ESD tracks
753         TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
754         Int_t nGoodTracks = list->GetEntries();
755         Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
756   
757         labelArr = new Int_t[nGoodTracks];
758         etaArr = new Float_t[nGoodTracks];
759         thirdDimArr = new Float_t[nGoodTracks];
760   
761         // loop over esd tracks
762         for (Int_t i=0; i<nGoodTracks; i++)
763         {
764           AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
765           if (!esdTrack)
766           {
767             AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
768             continue;
769           }
770           
771           Float_t phi = esdTrack->Phi();
772           if (phi < 0)
773             phi += TMath::Pi() * 2;
774   
775           Float_t eta = esdTrack->Eta();
776           Int_t label = TMath::Abs(esdTrack->GetLabel());
777           Float_t pT  = esdTrack->Pt();
778           
779           // force pT to fixed value without B field
780           if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
781             pT = 1;
782   
783           fPhi->Fill(phi);
784           fEtaPhi->Fill(eta, phi);
785           if (eventTriggered && vtxESD)
786             fRawPt->Fill(pT);
787   
788           if (esdTrack->Pt() < 0.15)
789             continue;
790           
791           // INEL>0 trigger
792           if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
793             foundInEta10 = kTRUE;
794   
795           if (fOnlyPrimaries)
796           {
797             if (label == 0)
798               continue;
799             
800             if (stack->IsPhysicalPrimary(label) == kFALSE)
801               continue;
802           }
803   
804           if (fUseMCKine)
805           {
806             if (label > 0)
807             {
808               TParticle* particle = stack->Particle(label);
809               eta = particle->Eta();
810               pT = particle->Pt();
811             }
812             else
813               Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
814           }
815   
816           if (fSymmetrize)
817             eta = TMath::Abs(eta);
818           etaArr[inputCount] = eta;
819           labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
820           thirdDimArr[inputCount] = pT;
821           ++inputCount;
822         }
823         
824         //if (inputCount > 30)
825         //  Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber());
826         
827         // TODO restrict inputCount used as measure for the multiplicity to |eta| < 1
828   
829         delete list;
830       }
831       
832       if (!foundInEta10)
833         eventTriggered = kFALSE;
834     }
835     else
836       return;
837
838     // Processing of ESD information (always)
839     if (eventTriggered)
840     {
841       // control hist
842       fMult->Fill(inputCount);
843       fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
844
845       if (vtxESD)
846       {
847         // control hist
848         
849         if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
850           fVertexVsMult->Fill(vtxESD->GetXv(), vtxESD->GetYv(), inputCount);
851       
852         Int_t part05 = 0;
853         Int_t part10 = 0;
854         for (Int_t i=0; i<inputCount; ++i)
855         {
856           Float_t eta = etaArr[i];
857           Float_t thirdDim = thirdDimArr[i];
858
859           fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, thirdDim);
860
861           if (TMath::Abs(vtx[2]) < 10)
862           {
863             fPartEta[0]->Fill(eta);
864
865             if (vtx[2] < 0)
866               fPartEta[1]->Fill(eta);
867             else
868               fPartEta[2]->Fill(eta);
869           }
870           
871           if (TMath::Abs(eta) < 0.5)
872             part05++;
873           if (TMath::Abs(eta) < 1.0)
874             part10++;
875         }
876         
877         Int_t multAxis = inputCount;
878         if (fMultAxisEta1)
879           multAxis = part10;
880
881         fMultVtx->Fill(multAxis);
882         //if (TMath::Abs(vtx[2]) < 10)
883         //  fMultVtx->Fill(part05);
884
885         // for event count per vertex
886         fdNdEtaAnalysisESD->FillEvent(vtx[2], multAxis);
887
888         // control hist
889         if (multAxis > 0)
890           fEvents->Fill(vtx[2]);
891
892         if (fReadMC)
893         {
894           // from tracks is only done for triggered and vertex reconstructed events
895           for (Int_t i=0; i<inputCount; ++i)
896           {
897             Int_t label = labelArr[i];
898
899             if (label < 0)
900               continue;
901
902             //Printf("Getting particle of track %d", label);
903             TParticle* particle = stack->Particle(label);
904
905             if (!particle)
906             {
907               AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
908               continue;
909             }
910
911             Float_t thirdDim = -1;
912             if (fAnalysisMode & AliPWG0Helper::kSPD)
913             {
914               if (fFillPhi)
915               {
916                 thirdDim = particle->Phi();
917               }
918               else
919                 thirdDim = multAxis;
920             }
921             else
922               thirdDim = particle->Pt();
923
924             Float_t eta = particle->Eta();
925             if (fSymmetrize)
926               eta = TMath::Abs(eta);
927             fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], eta, thirdDim);
928           } // end of track loop
929
930           // for event count per vertex
931           fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], multAxis);
932         }
933       }
934     }
935     if (etaArr)
936       delete[] etaArr;
937     if (labelArr)
938       delete[] labelArr;
939     if (thirdDimArr)
940       delete[] thirdDimArr;
941   }
942
943   if (fReadMC)   // Processing of MC information (optional)
944   {
945     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
946     if (!eventHandler) {
947       Printf("ERROR: Could not retrieve MC event handler");
948       return;
949     }
950
951     AliMCEvent* mcEvent = eventHandler->MCEvent();
952     if (!mcEvent) {
953       Printf("ERROR: Could not retrieve MC event");
954       return;
955     }
956
957     AliStack* stack = mcEvent->Stack();
958     if (!stack)
959     {
960       AliDebug(AliLog::kError, "Stack not available");
961       return;
962     }
963
964     AliHeader* header = mcEvent->Header();
965     if (!header)
966     {
967       AliDebug(AliLog::kError, "Header not available");
968       return;
969     }
970
971     // get the MC vertex
972     AliGenEventHeader* genHeader = header->GenEventHeader();
973     if (!genHeader)
974     {
975       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
976       return;
977     }
978
979     TArrayF vtxMC(3);
980     genHeader->PrimaryVertex(vtxMC);
981     
982     // get process type
983     Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
984     AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
985
986     if (processType == AliPWG0Helper::kInvalidProcess)
987       AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
988
989     // loop over mc particles
990     Int_t nPrim  = stack->GetNprimary();
991
992     Int_t nAcceptedParticles = 0;
993     Bool_t oneParticleEvent = kFALSE;
994
995     // count particles first, then fill
996     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
997     {
998       if (!stack->IsPhysicalPrimary(iMc))
999         continue;
1000
1001       //Printf("Getting particle %d", iMc);
1002       TParticle* particle = stack->Particle(iMc);
1003
1004       if (!particle)
1005         continue;
1006
1007       if (particle->GetPDG()->Charge() == 0)
1008         continue;
1009
1010       //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
1011       Float_t eta = particle->Eta();
1012
1013       if (TMath::Abs(eta) < 1.0)
1014         oneParticleEvent = kTRUE;
1015
1016       // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: these histograms are only gathered for comparison))
1017       if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
1018         nAcceptedParticles++;
1019     }
1020
1021     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
1022     {
1023       if (!stack->IsPhysicalPrimary(iMc))
1024         continue;
1025
1026       //Printf("Getting particle %d", iMc);
1027       TParticle* particle = stack->Particle(iMc);
1028
1029       if (!particle)
1030         continue;
1031
1032       if (particle->GetPDG()->Charge() == 0)
1033         continue;
1034
1035       Float_t eta = particle->Eta();
1036       if (fSymmetrize)
1037         eta = TMath::Abs(eta);
1038
1039       Float_t thirdDim = -1;
1040
1041       if (fAnalysisMode & AliPWG0Helper::kSPD)
1042       {
1043         if (fFillPhi)
1044         {
1045           thirdDim = particle->Phi();
1046         }
1047         else
1048           thirdDim = nAcceptedParticles;
1049       }
1050       else
1051         thirdDim = particle->Pt();
1052
1053       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
1054
1055       if (processType != AliPWG0Helper::kSD)
1056         fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
1057
1058       if (processType == AliPWG0Helper::kND)
1059         fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
1060       
1061       if (oneParticleEvent)
1062         fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
1063
1064       if (eventTriggered)
1065       {
1066         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
1067         if (vtxESD)
1068           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
1069       }
1070
1071       if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0 && particle->P() > 0)
1072       {
1073         //Float_t value = 1. / TMath::TwoPi() / particle->Pt() * particle->Energy() / particle->P();
1074         Float_t value = 1;
1075         fPartPt->Fill(particle->Pt(), value);
1076       }
1077     }
1078
1079     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
1080     if (processType != AliPWG0Helper::kSD)
1081       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
1082     if (processType == AliPWG0Helper::kND)
1083       fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
1084     if (oneParticleEvent)
1085       fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles);
1086
1087     if (eventTriggered)
1088     {
1089       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
1090       if (vtxESD)
1091         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
1092     }
1093   }
1094 }
1095
1096 void AlidNdEtaTask::Terminate(Option_t *)
1097 {
1098   // The Terminate() function is the last function to be called during
1099   // a query. It always runs on the client, it can be used to present
1100   // the results graphically or save the results to file.
1101
1102   fOutput = dynamic_cast<TList*> (GetOutputData(0));
1103   if (!fOutput)
1104     Printf("ERROR: fOutput not available");
1105
1106   if (fOutput)
1107   {
1108     fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
1109     fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
1110     fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
1111     for (Int_t i=0; i<3; ++i)
1112       fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
1113     fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
1114     fVertexResolution = dynamic_cast<TH2F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
1115
1116     fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
1117     fVertexVsMult = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexVsMult"));
1118     fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
1119     fRawPt = dynamic_cast<TH1F*> (fOutput->FindObject("fRawPt"));
1120     fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
1121     for (Int_t i=0; i<2; ++i)
1122       fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
1123     fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
1124     fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
1125     fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
1126     fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
1127     fTrackletsVsClusters = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsClusters"));
1128     fTrackletsVsUnassigned = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsUnassigned"));
1129     fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
1130     fStats2 = dynamic_cast<TH2F*> (fOutput->FindObject("fStats2"));
1131
1132     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
1133   }
1134
1135   if (!fdNdEtaAnalysisESD)
1136   {
1137     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
1138   }
1139   else
1140   {
1141     if (fMult && fMultVtx)
1142     {
1143       new TCanvas;
1144       fMult->Draw();
1145       fMultVtx->SetLineColor(2);
1146       fMultVtx->Draw("SAME");
1147     }
1148
1149     if (fFiredChips)
1150     {
1151       new TCanvas;
1152       fFiredChips->Draw("COLZ");
1153     }
1154
1155     if (fPartEta[0])
1156     {
1157       Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-9.9), fEvents->GetXaxis()->FindBin(-0.001));
1158       Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(9.9));
1159
1160       Printf("%d events with vertex used", events1 + events2);
1161       Printf("%d total events with vertex, %d total triggered events, %d triggered events with 0 multiplicity", (Int_t) fMultVtx->Integral(), (Int_t) fMult->Integral(), (Int_t) fMult->GetBinContent(1));
1162
1163       if (events1 > 0 && events2 > 0)
1164       {
1165         fPartEta[0]->Scale(1.0 / (events1 + events2));
1166         fPartEta[1]->Scale(1.0 / events1);
1167         fPartEta[2]->Scale(1.0 / events2);
1168
1169         for (Int_t i=0; i<3; ++i)
1170           fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
1171
1172         new TCanvas("control", "control", 500, 500);
1173         for (Int_t i=0; i<3; ++i)
1174         {
1175           fPartEta[i]->SetLineColor(i+1);
1176           fPartEta[i]->Draw((i==0) ? "" : "SAME");
1177         }
1178       }
1179     }
1180
1181     if (fEvents)
1182     {
1183         new TCanvas("control3", "control3", 500, 500);
1184         fEvents->Draw();
1185     }
1186     
1187     if (fStats2)
1188     {
1189       new TCanvas;
1190       fStats2->Draw("TEXT");
1191     }
1192
1193     TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
1194
1195     if (fdNdEtaAnalysisESD)
1196       fdNdEtaAnalysisESD->SaveHistograms();
1197
1198     if (fEsdTrackCuts)
1199       fEsdTrackCuts->SaveHistograms("esd_track_cuts");
1200
1201     if (fMult)
1202       fMult->Write();
1203
1204     if (fMultVtx)
1205       fMultVtx->Write();
1206
1207     for (Int_t i=0; i<3; ++i)
1208       if (fPartEta[i])
1209         fPartEta[i]->Write();
1210
1211     if (fEvents)
1212       fEvents->Write();
1213
1214     if (fVertexResolution)
1215     {
1216       fVertexResolution->Write();
1217       fVertexResolution->ProjectionX();
1218       fVertexResolution->ProjectionY();
1219     }
1220
1221     if (fDeltaPhi)
1222       fDeltaPhi->Write();
1223
1224     if (fDeltaTheta)
1225       fDeltaTheta->Write();
1226     
1227     if (fPhi)
1228       fPhi->Write();
1229
1230     if (fRawPt)
1231       fRawPt->Write();
1232
1233     if (fEtaPhi)
1234       fEtaPhi->Write();
1235
1236     for (Int_t i=0; i<2; ++i)
1237       if (fZPhi[i])
1238         fZPhi[i]->Write();
1239     
1240     if (fModuleMap)
1241       fModuleMap->Write();
1242     
1243     if (fFiredChips)
1244       fFiredChips->Write();
1245
1246     if (fTrackletsVsClusters)
1247       fTrackletsVsClusters->Write();
1248     
1249     if (fTrackletsVsUnassigned)
1250       fTrackletsVsUnassigned->Write();
1251     
1252     if (fStats)
1253       fStats->Write();
1254
1255     if (fStats2)
1256       fStats2->Write();
1257     
1258     if (fVertex)
1259       fVertex->Write();
1260
1261     if (fVertexVsMult)
1262       fVertexVsMult->Write();
1263     
1264     fout->Write();
1265     fout->Close();
1266
1267     Printf("Writting result to analysis_esd_raw.root");
1268   }
1269
1270   if (fReadMC)
1271   {
1272     if (fOutput)
1273     {
1274       fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
1275       fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
1276       fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
1277       fdNdEtaAnalysisOnePart = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaOnePart"));
1278       fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
1279       fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
1280       fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
1281       fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
1282     }
1283
1284     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
1285     {
1286       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
1287       return;
1288     }
1289
1290     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
1291     fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
1292     fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
1293     fdNdEtaAnalysisOnePart->Finish(0, -1, AlidNdEtaCorrection::kNone);
1294     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
1295     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
1296     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
1297
1298     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
1299     fPartPt->Scale(1.0/events);
1300     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
1301
1302     TFile* fout = new TFile("analysis_mc.root","RECREATE");
1303
1304     fdNdEtaAnalysis->SaveHistograms();
1305     fdNdEtaAnalysisND->SaveHistograms();
1306     fdNdEtaAnalysisNSD->SaveHistograms();
1307     fdNdEtaAnalysisOnePart->SaveHistograms();
1308     fdNdEtaAnalysisTr->SaveHistograms();
1309     fdNdEtaAnalysisTrVtx->SaveHistograms();
1310     fdNdEtaAnalysisTracks->SaveHistograms();
1311
1312     if (fPartPt)
1313       fPartPt->Write();
1314
1315     fout->Write();
1316     fout->Close();
1317
1318     Printf("Writting result to analysis_mc.root");
1319
1320     if (fPartPt)
1321     {
1322       new TCanvas("control2", "control2", 500, 500);
1323       fPartPt->Draw();
1324     }
1325   }
1326 }