66f0d09d9bb835da717036e431464667d26d48b5
[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
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 "AliCorrection.h"
35 #include "AliCorrectionMatrix3D.h"
36 #include "dNdEta/dNdEtaAnalysis.h"
37
38 ClassImp(AlidNdEtaTask)
39
40 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
41   AliAnalysisTask("AlidNdEtaTask", ""),
42   fESD(0),
43   fOutput(0),
44   fOption(opt),
45   fAnalysisMode(AliPWG0Helper::kTPC),
46   fTrigger(AliPWG0Helper::kMB1),
47   fReadMC(kFALSE),
48   fUseMCVertex(kFALSE),
49   fOnlyPrimaries(kFALSE),
50   fUseMCKine(kFALSE),
51   fEsdTrackCuts(0),
52   fdNdEtaAnalysisESD(0),
53   fMult(0),
54   fMultVtx(0),
55   fEvents(0),
56   fVertexResolution(0),
57   fdNdEtaAnalysis(0),
58   fdNdEtaAnalysisNSD(0),
59   fdNdEtaAnalysisTr(0),
60   fdNdEtaAnalysisTrVtx(0),
61   fdNdEtaAnalysisTracks(0),
62   fVertex(0),
63   fPartPt(0),
64   fPhi(0),
65   fEtaPhi(0),
66   fDeltaPhi(0)
67 {
68   //
69   // Constructor. Initialization of pointers
70   //
71
72   // Define input and output slots here
73   DefineInput(0, TChain::Class());
74   DefineOutput(0, TList::Class());
75
76   AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
77 }
78
79 AlidNdEtaTask::~AlidNdEtaTask()
80 {
81   //
82   // Destructor
83   //
84
85   // histograms are in the output list and deleted when the output
86   // list is deleted by the TSelector dtor
87
88   if (fOutput) {
89     delete fOutput;
90     fOutput = 0;
91   }
92 }
93
94 Bool_t AlidNdEtaTask::Notify()
95 {
96   static Int_t count = 0;
97   count++;
98   Printf("Processing %d. file", count);
99   return kTRUE;
100 }
101
102 //________________________________________________________________________
103 void AlidNdEtaTask::ConnectInputData(Option_t *)
104 {
105   // Connect ESD
106   // Called once
107
108   Printf("AlidNdEtaTask::ConnectInputData called");
109
110   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
111
112   if (!esdH) {
113     Printf("ERROR: Could not get ESDInputHandler");
114   } else {
115     fESD = esdH->GetEvent();
116
117     // Enable only the needed branches
118     esdH->SetActiveBranches("AliESDHeader Vertex");
119
120     if (fAnalysisMode == AliPWG0Helper::kSPD)
121       esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
122
123     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
124       esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
125     }
126   }
127
128   // disable info messages of AliMCEvent (per event)
129   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
130 }
131
132 void AlidNdEtaTask::CreateOutputObjects()
133 {
134   // create result objects and add to output list
135
136   Printf("AlidNdEtaTask::CreateOutputObjects");
137
138   fOutput = new TList;
139   fOutput->SetOwner();
140
141   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
142   fOutput->Add(fdNdEtaAnalysisESD);
143
144   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
145   fOutput->Add(fMult);
146
147   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
148   fOutput->Add(fMultVtx);
149
150   for (Int_t i=0; i<3; ++i)
151   {
152     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
153     fPartEta[i]->Sumw2();
154     fOutput->Add(fPartEta[i]);
155   }
156
157   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 160, -40, 40);
158   fOutput->Add(fEvents);
159
160   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
161   fOutput->Add(fVertexResolution);
162
163   fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
164   fOutput->Add(fPhi);
165
166   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*5, 0, 2 * TMath::Pi());
167   fOutput->Add(fEtaPhi);
168
169   if (fAnalysisMode == AliPWG0Helper::kSPD)
170     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 18*50, -3.14, 3.14);
171
172   if (fReadMC)
173   {
174     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
175     fOutput->Add(fdNdEtaAnalysis);
176
177     fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
178     fOutput->Add(fdNdEtaAnalysisNSD);
179
180     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
181     fOutput->Add(fdNdEtaAnalysisTr);
182
183     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
184     fOutput->Add(fdNdEtaAnalysisTrVtx);
185
186     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
187     fOutput->Add(fdNdEtaAnalysisTracks);
188
189     fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
190     fOutput->Add(fVertex);
191
192     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
193     fPartPt->Sumw2();
194     fOutput->Add(fPartPt);
195   }
196
197   if (fEsdTrackCuts)
198   {
199     fEsdTrackCuts->SetName("fEsdTrackCuts");
200     fOutput->Add(fEsdTrackCuts);
201   }
202 }
203
204 void AlidNdEtaTask::Exec(Option_t*)
205 {
206   // process the event
207
208   // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
209   Bool_t eventTriggered = kFALSE;
210   const AliESDVertex* vtxESD = 0;
211
212   // ESD analysis
213   if (fESD)
214   {
215     // trigger definition
216     eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);
217
218     // get the ESD vertex
219     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
220
221     Double_t vtx[3];
222     if (vtxESD) {
223       vtxESD->GetXYZ(vtx);
224     }
225     else
226       Printf("No vertex");
227
228     // fill z vertex resolution
229     if (vtxESD)
230       fVertexResolution->Fill(vtxESD->GetZRes());
231
232     // post the data already here
233     PostData(0, fOutput);
234
235     // needed for syst. studies
236     AliStack* stack = 0;
237     TArrayF vtxMC(3);
238
239     if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
240       if (!fReadMC) {
241         Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
242         return;
243       }
244
245       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
246       if (!eventHandler) {
247         Printf("ERROR: Could not retrieve MC event handler");
248         return;
249       }
250
251       AliMCEvent* mcEvent = eventHandler->MCEvent();
252       if (!mcEvent) {
253         Printf("ERROR: Could not retrieve MC event");
254         return;
255       }
256
257       AliHeader* header = mcEvent->Header();
258       if (!header)
259       {
260         AliDebug(AliLog::kError, "Header not available");
261         return;
262       }
263
264       // get the MC vertex
265       AliGenEventHeader* genHeader = header->GenEventHeader();
266       if (!genHeader)
267       {
268         AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
269         return;
270       }
271       genHeader->PrimaryVertex(vtxMC);
272
273       if (fUseMCVertex)
274       {
275         Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
276         vtx[2] = vtxMC[2];
277       }
278
279       stack = mcEvent->Stack();
280       if (!stack)
281       {
282         AliDebug(AliLog::kError, "Stack not available");
283         return;
284       }
285     }
286
287     // create list of (label, eta, pt) tuples
288     Int_t inputCount = 0;
289     Int_t* labelArr = 0;
290     Float_t* etaArr = 0;
291     Float_t* ptArr = 0;
292     if (fAnalysisMode == AliPWG0Helper::kSPD)
293     {
294       // get tracklets
295       const AliMultiplicity* mult = fESD->GetMultiplicity();
296       if (!mult)
297       {
298         AliDebug(AliLog::kError, "AliMultiplicity not available");
299         return;
300       }
301
302       labelArr = new Int_t[mult->GetNumberOfTracklets()];
303       etaArr = new Float_t[mult->GetNumberOfTracklets()];
304       ptArr = new Float_t[mult->GetNumberOfTracklets()];
305
306       if (fOnlyPrimaries)
307         Printf("Processing only primaries (MC information used). This is for systematical checks only.");
308
309       // get multiplicity from ITS tracklets
310       for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
311       {
312         //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
313
314         if (fOnlyPrimaries)
315           if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
316             continue;
317
318         Float_t deltaPhi = mult->GetDeltaPhi(i);
319         // prevent values to be shifted by 2 Pi()
320         if (deltaPhi < -TMath::Pi())
321           deltaPhi += TMath::Pi() * 2;
322         if (deltaPhi > TMath::Pi())
323           deltaPhi -= TMath::Pi() * 2;
324
325         if (TMath::Abs(deltaPhi) > 1)
326           printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
327
328         Float_t phi = mult->GetPhi(i);
329         if (phi < 0)
330           phi += TMath::Pi() * 2;
331         fPhi->Fill(phi);
332         fEtaPhi->Fill(mult->GetEta(i), phi);
333
334         fDeltaPhi->Fill(deltaPhi);
335
336         // TODO implement fUseMCKine here
337
338         etaArr[inputCount] = mult->GetEta(i);
339         labelArr[inputCount] = mult->GetLabel(i, 0);
340         ptArr[inputCount] = 0; // no pt for tracklets
341         ++inputCount;
342       }
343
344       Printf("Accepted %d tracklets", inputCount);
345     }
346     else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
347     {
348       if (!fEsdTrackCuts)
349       {
350         AliDebug(AliLog::kError, "fESDTrackCuts not available");
351         return;
352       }
353
354       // get multiplicity from ESD tracks
355       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
356       Int_t nGoodTracks = list->GetEntries();
357
358       labelArr = new Int_t[nGoodTracks];
359       etaArr = new Float_t[nGoodTracks];
360       ptArr = new Float_t[nGoodTracks];
361
362       if (fOnlyPrimaries)
363         Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
364
365       if (fUseMCKine)
366         Printf("WARNING: Using MC kine information. This is for systematical checks only.");
367
368       // loop over esd tracks
369       for (Int_t i=0; i<nGoodTracks; i++)
370       {
371         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
372         if (!esdTrack)
373         {
374           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
375           continue;
376         }
377
378         Float_t phi = esdTrack->Phi();
379         if (phi < 0)
380           phi += TMath::Pi() * 2;
381         fPhi->Fill(phi);
382         fEtaPhi->Fill(esdTrack->Eta(), phi);
383
384         Float_t eta = esdTrack->Eta();
385         Int_t label = TMath::Abs(esdTrack->GetLabel());
386         Float_t pT  = esdTrack->Pt();
387
388         if (fOnlyPrimaries && label == 0)
389           continue;
390
391         if (fUseMCKine)
392         {
393           if (label > 0)
394           {
395             TParticle* particle = stack->Particle(label);
396             eta = particle->Eta();
397             pT = particle->Pt();
398           }
399           else
400             Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
401         }
402
403         etaArr[inputCount] = eta;
404         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
405         ptArr[inputCount] = pT;
406         ++inputCount;
407       }
408
409       Printf("Accepted %d tracks", nGoodTracks);
410
411       delete list;
412     }
413     else
414       return;
415
416     // Processing of ESD information (always)
417     if (eventTriggered)
418     {
419       // control hist
420       fMult->Fill(inputCount);
421       fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
422
423       if (vtxESD)
424       {
425         // control hist
426         fMultVtx->Fill(inputCount);
427
428         for (Int_t i=0; i<inputCount; ++i)
429         {
430           Float_t eta = etaArr[i];
431           Float_t pt  = ptArr[i];
432
433           fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
434
435           if (TMath::Abs(vtx[2]) < 20)
436           {
437             fPartEta[0]->Fill(eta);
438
439             if (vtx[2] < 0)
440               fPartEta[1]->Fill(eta);
441             else
442               fPartEta[2]->Fill(eta);
443           }
444         }
445
446         // for event count per vertex
447         fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
448
449         // control hist
450         fEvents->Fill(vtx[2]);
451
452         if (fReadMC)
453         {
454           // from tracks is only done for triggered and vertex reconstructed events
455           for (Int_t i=0; i<inputCount; ++i)
456           {
457             Int_t label = labelArr[i];
458
459             if (label < 0)
460               continue;
461
462             //Printf("Getting particle of track %d", label);
463             TParticle* particle = stack->Particle(label);
464
465             if (!particle)
466             {
467               AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
468               continue;
469             }
470
471             fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
472           } // end of track loop
473
474           // for event count per vertex
475           fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
476         }
477       }
478     }
479
480     delete[] etaArr;
481     delete[] labelArr;
482     delete[] ptArr;
483   }
484
485   if (fReadMC)   // Processing of MC information (optional)
486   {
487     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
488     if (!eventHandler) {
489       Printf("ERROR: Could not retrieve MC event handler");
490       return;
491     }
492
493     AliMCEvent* mcEvent = eventHandler->MCEvent();
494     if (!mcEvent) {
495       Printf("ERROR: Could not retrieve MC event");
496       return;
497     }
498
499     AliStack* stack = mcEvent->Stack();
500     if (!stack)
501     {
502       AliDebug(AliLog::kError, "Stack not available");
503       return;
504     }
505
506     AliHeader* header = mcEvent->Header();
507     if (!header)
508     {
509       AliDebug(AliLog::kError, "Header not available");
510       return;
511     }
512
513     // get the MC vertex
514     AliGenEventHeader* genHeader = header->GenEventHeader();
515     if (!genHeader)
516     {
517       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
518       return;
519     }
520
521     TArrayF vtxMC(3);
522     genHeader->PrimaryVertex(vtxMC);
523
524     // get process type
525     Int_t processType = AliPWG0Helper::GetEventProcessType(header);
526     AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
527
528     if (processType==AliPWG0Helper::kInvalidProcess)
529       AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
530
531     // loop over mc particles
532     Int_t nPrim  = stack->GetNprimary();
533
534     Int_t nAcceptedParticles = 0;
535
536     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
537     {
538       if (!stack->IsPhysicalPrimary(iMc))
539         continue;
540
541       //Printf("Getting particle %d", iMc);
542       TParticle* particle = stack->Particle(iMc);
543
544       if (!particle)
545         continue;
546
547       if (particle->GetPDG()->Charge() == 0)
548         continue;
549
550       //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
551       Float_t eta = particle->Eta();
552       Float_t pt = particle->Pt();
553
554        // 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))
555       if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
556         nAcceptedParticles++;
557
558       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
559       fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
560
561       if (processType != AliPWG0Helper::kSD )
562         fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, pt);
563
564       if (eventTriggered)
565       {
566         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
567         if (vtxESD)
568           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
569       }
570
571       if (TMath::Abs(eta) < 0.8)
572         fPartPt->Fill(pt);
573     }
574
575     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
576     if (processType != AliPWG0Helper::kSD)
577       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
578
579     if (eventTriggered)
580     {
581       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
582       if (vtxESD)
583         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
584     }
585   }
586 }
587
588 void AlidNdEtaTask::Terminate(Option_t *)
589 {
590   // The Terminate() function is the last function to be called during
591   // a query. It always runs on the client, it can be used to present
592   // the results graphically or save the results to file.
593
594   fOutput = dynamic_cast<TList*> (GetOutputData(0));
595   if (!fOutput)
596     Printf("ERROR: fOutput not available");
597
598   if (fOutput)
599   {
600     fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
601     fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
602     fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
603     for (Int_t i=0; i<3; ++i)
604       fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
605     fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
606     fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
607
608     fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
609     fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
610     fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
611
612     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
613   }
614
615   if (!fdNdEtaAnalysisESD)
616   {
617     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
618     return;
619   }
620
621   if (fMult && fMultVtx)
622   {
623     new TCanvas;
624     fMult->Draw();
625     fMultVtx->SetLineColor(2);
626     fMultVtx->Draw("SAME");
627   }
628
629   if (fPartEta[0])
630   {
631     Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
632     Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
633
634     Printf("%d events with vertex used", events1 + events2);
635
636     if (events1 > 0 && events2 > 0)
637     {
638     fPartEta[0]->Scale(1.0 / (events1 + events2));
639     fPartEta[1]->Scale(1.0 / events1);
640     fPartEta[2]->Scale(1.0 / events2);
641
642     for (Int_t i=0; i<3; ++i)
643       fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
644
645     new TCanvas("control", "control", 500, 500);
646     for (Int_t i=0; i<3; ++i)
647     {
648       fPartEta[i]->SetLineColor(i+1);
649       fPartEta[i]->Draw((i==0) ? "" : "SAME");
650     }
651     }
652   }
653
654   if (fEvents)
655   {
656     new TCanvas("control3", "control3", 500, 500);
657     fEvents->Draw();
658   }
659
660   TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
661
662   if (fdNdEtaAnalysisESD)
663     fdNdEtaAnalysisESD->SaveHistograms();
664
665   if (fEsdTrackCuts)
666     fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
667
668   if (fMult)
669     fMult->Write();
670
671   if (fMultVtx)
672     fMultVtx->Write();
673
674   for (Int_t i=0; i<3; ++i)
675     if (fPartEta[i])
676       fPartEta[i]->Write();
677
678   if (fEvents)
679     fEvents->Write();
680
681   if (fVertexResolution)
682     fVertexResolution->Write();
683
684   if (fDeltaPhi)
685     fDeltaPhi->Write();
686
687   if (fPhi)
688     fPhi->Write();
689
690   if (fEtaPhi)
691     fEtaPhi->Write();
692
693   fout->Write();
694   fout->Close();
695
696   Printf("Writting result to analysis_esd_raw.root");
697
698   if (fReadMC)
699   {
700     if (fOutput)
701     {
702       fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
703       fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
704       fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
705       fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
706       fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
707       fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
708       fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
709     }
710
711     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
712     {
713       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
714       return;
715     }
716
717     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
718     fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
719     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
720     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
721     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
722
723     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
724     fPartPt->Scale(1.0/events);
725     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
726
727     fout = new TFile("analysis_mc.root","RECREATE");
728
729     fdNdEtaAnalysis->SaveHistograms();
730     fdNdEtaAnalysisNSD->SaveHistograms();
731     fdNdEtaAnalysisTr->SaveHistograms();
732     fdNdEtaAnalysisTrVtx->SaveHistograms();
733     fdNdEtaAnalysisTracks->SaveHistograms();
734
735     if (fPartPt)
736       fPartPt->Write();
737
738     if (fVertex)
739       fVertex->Write();
740
741     fout->Write();
742     fout->Close();
743
744     Printf("Writting result to analysis_mc.root");
745
746     if (fPartPt)
747     {
748       new TCanvas("control2", "control2", 500, 500);
749       fPartPt->Draw();
750     }
751   }
752 }