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