0bc47389e86eb8c2d4d8723a16f932c84d2b1a4c
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMCCorrectionsTask.cxx
1 // 
2 // Calculate the corrections in the forward regions
3 // 
4 // Inputs: 
5 //   - AliESDEvent 
6 //
7 // Outputs: 
8 //   - AliAODForwardMult 
9 // 
10 // Histograms 
11 //   
12 // Corrections used 
13 // 
14 #include "AliForwardMCCorrectionsTask.h"
15 #include "AliTriggerAnalysis.h"
16 #include "AliPhysicsSelection.h"
17 #include "AliLog.h"
18 #include "AliHeader.h"
19 #include "AliGenEventHeader.h"
20 #include "AliESDEvent.h"
21 #include "AliAODHandler.h"
22 #include "AliMultiplicity.h"
23 #include "AliInputEventHandler.h"
24 #include "AliStack.h"
25 #include "AliMCEvent.h"
26 #include "AliAODForwardMult.h"
27 #include "AliFMDStripIndex.h"
28 #include "AliFMDCorrSecondaryMap.h"
29 #include <TH1.h>
30 #include <TH2D.h>
31 #include <TDirectory.h>
32 #include <TTree.h>
33 #include <TList.h>
34 #include <TROOT.h>
35 #include <iostream>
36
37 //====================================================================
38 namespace {
39   const char* GetEventName(Bool_t tr, Bool_t vtx) 
40   {
41     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
42   }
43   const char* GetHitsName(UShort_t d, Char_t r) 
44   {
45     return Form("hitsFMD%d%c", d, r);
46   }
47   const char* GetStripsName(UShort_t d, Char_t r)
48   {
49     return Form("stripsFMD%d%c", d, r);
50   }
51   const char* GetPrimaryName(Char_t r, Bool_t trVtx)
52   {
53     return Form("primaries%s%s", 
54                 ((r == 'I' || r == 'i') ? "Inner" : "Outer"), 
55                 (trVtx ? "TrVtx" : "All"));
56   }
57 }
58
59 //====================================================================
60 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
61   : AliAnalysisTaskSE(),
62     fInspector(),
63     fTrackDensity(),
64     fESDFMD(),
65     fVtxBins(0),
66     fFirstEvent(true),
67     fHEvents(0), 
68     fHEventsTr(0), 
69     fHEventsTrVtx(0),
70     fVtxAxis(),
71     fEtaAxis(),
72     fList()
73 {
74   // 
75   // Constructor 
76   // 
77   // Parameters:
78   //    name Name of task 
79   //
80 }
81
82 //____________________________________________________________________
83 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
84   : AliAnalysisTaskSE(name),
85     fInspector("eventInspector"), 
86     fTrackDensity("trackDensity"),
87     fESDFMD(),
88     fVtxBins(0),
89     fFirstEvent(true),
90     fHEvents(0), 
91     fHEventsTr(0), 
92     fHEventsTrVtx(0),
93     fVtxAxis(10,-10,10), 
94     fEtaAxis(200,-4,6),
95     fList()
96 {
97   // 
98   // Constructor 
99   // 
100   // Parameters:
101   //    name Name of task 
102   //
103   DefineOutput(1, TList::Class());
104   DefineOutput(2, TList::Class());
105 }
106
107 //____________________________________________________________________
108 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
109   : AliAnalysisTaskSE(o),
110     fInspector(o.fInspector),
111     fTrackDensity(),
112     fESDFMD(o.fESDFMD),
113     fVtxBins(0),
114     fFirstEvent(o.fFirstEvent),
115     fHEvents(o.fHEvents), 
116     fHEventsTr(o.fHEventsTr), 
117     fHEventsTrVtx(o.fHEventsTrVtx),
118     fVtxAxis(10,-10,10), 
119     fEtaAxis(200,-4,6),
120     fList(o.fList)
121 {
122   // 
123   // Copy constructor 
124   // 
125   // Parameters:
126   //    o Object to copy from 
127   //
128 }
129
130 //____________________________________________________________________
131 AliForwardMCCorrectionsTask&
132 AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
133 {
134   // 
135   // Assignment operator 
136   // 
137   // Parameters:
138   //    o Object to assign from 
139   // 
140   // Return:
141   //    Reference to this object 
142   //
143   fInspector         = o.fInspector;
144   fTrackDensity      = o.fTrackDensity;
145   fESDFMD            = o.fESDFMD;
146   fVtxBins           = o.fVtxBins;
147   fFirstEvent        = o.fFirstEvent;
148   fHEvents           = o.fHEvents;
149   fHEventsTr         = o.fHEventsTr;
150   fHEventsTrVtx      = o.fHEventsTrVtx;
151   SetVertexAxis(o.fVtxAxis);
152   SetEtaAxis(o.fEtaAxis);
153
154   return *this;
155 }
156
157 //____________________________________________________________________
158 void
159 AliForwardMCCorrectionsTask::Init()
160 {
161   // 
162   // Initialize the task 
163   // 
164   //
165 }
166
167 //____________________________________________________________________
168 void
169 AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min, 
170                                            Double_t max)
171 {
172   // 
173   // Set the vertex axis to use
174   // 
175   // Parameters:
176   //    nBins Number of bins
177   //    vzMin Least @f$z@f$ coordinate of interation point
178   //    vzMax Largest @f$z@f$ coordinate of interation point
179   //
180   if (max < min) max = -min;
181   if (min < max) { 
182     Double_t tmp = min;
183     min          = max;
184     max          = tmp;
185   }
186   if (min < -10) 
187     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
188   if (max > +10) 
189     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
190   fVtxAxis.Set(nBin, min, max);
191 }
192 //____________________________________________________________________
193 void
194 AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
195 {
196   // 
197   // Set the vertex axis to use
198   // 
199   // Parameters:
200   //    axis Axis
201   //
202   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
203 }
204
205 //____________________________________________________________________
206 void
207 AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
208 {
209   // 
210   // Set the eta axis to use
211   // 
212   // Parameters:
213   //    nBins Number of bins
214   //    vzMin Least @f$\eta@f$ 
215   //    vzMax Largest @f$\eta@f$ 
216   //
217   if (max < min) max = -min;
218   if (min < max) { 
219     Double_t tmp = min;
220     min          = max;
221     max          = tmp;
222   }
223   if (min < -4) 
224     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
225   if (max > +6) 
226     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
227   fVtxAxis.Set(nBin, min, max);
228 }
229 //____________________________________________________________________
230 void
231 AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
232 {
233   // 
234   // Set the eta axis to use
235   // 
236   // Parameters:
237   //    axis Axis
238   //
239   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
240 }
241
242 //____________________________________________________________________
243 void
244 AliForwardMCCorrectionsTask::DefineBins(TList* l)
245 {
246   if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
247   if (fVtxBins->GetEntries() > 0) return;
248
249   fVtxBins->SetOwner();
250   for (Int_t  i    = 1; i <= fVtxAxis.GetNbins(); i++) { 
251     Double_t  low  = fVtxAxis.GetBinLowEdge(i);
252     Double_t  high = fVtxAxis.GetBinUpEdge(i);
253     VtxBin*   bin  = new VtxBin(low, high, fEtaAxis);
254     fVtxBins->AddAt(bin, i);
255     bin->DefineOutput(l);
256   }
257 }
258
259 //____________________________________________________________________
260 void
261 AliForwardMCCorrectionsTask::UserCreateOutputObjects()
262 {
263   // 
264   // Create output objects 
265   // 
266   //
267   fList = new TList;
268   fList->SetOwner();
269   fList->SetName(Form("%sSums", GetName()));
270
271   DefineBins(fList);
272
273   fHEvents = new TH1I(GetEventName(false,false),
274                       "Number of all events", 
275                       fVtxAxis.GetNbins(), 
276                       fVtxAxis.GetXmin(), 
277                       fVtxAxis.GetXmax());
278   fHEvents->SetXTitle("v_{z} [cm]");
279   fHEvents->SetYTitle("# of events");
280   fHEvents->SetFillColor(kBlue+1);
281   fHEvents->SetFillStyle(3001);
282   fHEvents->SetDirectory(0);
283   fList->Add(fHEvents);
284
285   fHEventsTr = new TH1I(GetEventName(true, false), 
286                         "Number of triggered events",
287                         fVtxAxis.GetNbins(), 
288                         fVtxAxis.GetXmin(), 
289                         fVtxAxis.GetXmax());
290   fHEventsTr->SetXTitle("v_{z} [cm]");
291   fHEventsTr->SetYTitle("# of events");
292   fHEventsTr->SetFillColor(kRed+1);
293   fHEventsTr->SetFillStyle(3001);
294   fHEventsTr->SetDirectory(0);
295   fList->Add(fHEventsTr);
296
297   fHEventsTrVtx = new TH1I(GetEventName(true,true),
298                            "Number of events w/trigger and vertex", 
299                            fVtxAxis.GetNbins(), 
300                            fVtxAxis.GetXmin(), 
301                            fVtxAxis.GetXmax());
302   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
303   fHEventsTrVtx->SetYTitle("# of events");
304   fHEventsTrVtx->SetFillColor(kBlue+1);
305   fHEventsTrVtx->SetFillStyle(3001);
306   fHEventsTrVtx->SetDirectory(0);
307   fList->Add(fHEventsTrVtx);
308
309   // Copy axis objects to output 
310   TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis", 
311                           fVtxAxis.GetNbins(), 
312                           fVtxAxis.GetXmin(), 
313                           fVtxAxis.GetXmax());
314   TH1* etaAxis = new TH1D("etaAxis", "Eta axis", 
315                           fEtaAxis.GetNbins(), 
316                           fEtaAxis.GetXmin(), 
317                           fEtaAxis.GetXmax());
318   fList->Add(vtxAxis);
319   fList->Add(etaAxis);
320
321   AliInfo(Form("Initialising sub-routines: %p, %p", 
322                &fInspector, &fTrackDensity));
323   fInspector.DefineOutput(fList);
324   fInspector.Init(fVtxAxis);
325   fTrackDensity.DefineOutput(fList);
326
327   PostData(1, fList);
328 }
329 //____________________________________________________________________
330 void
331 AliForwardMCCorrectionsTask::UserExec(Option_t*)
332 {
333   // 
334   // Process each event 
335   // 
336   // Parameters:
337   //    option Not used
338   //  
339
340   // Get the input data - MC event
341   AliMCEvent*  mcEvent = MCEvent();
342   if (!mcEvent) { 
343     AliWarning("No MC event found");
344     return;
345   }
346
347   // Get the input data - ESD event
348   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
349   if (!esd) { 
350     AliWarning("No ESD event found for input event");
351     return;
352   }
353
354   //--- Read run information -----------------------------------------
355   if (fFirstEvent && esd->GetESDRun()) {
356     fInspector.ReadRunDetails(esd);
357     
358     AliInfo(Form("Initializing with parameters from the ESD:\n"
359                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
360                  "         AliESDEvent::GetBeamType()     ->%s\n"
361                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
362                  "         AliESDEvent::GetMagneticField()->%f\n"
363                  "         AliESDEvent::GetRunNumber()    ->%d\n",
364                  esd->GetBeamEnergy(),
365                  esd->GetBeamType(),
366                  esd->GetCurrentL3(),
367                  esd->GetMagneticField(),
368                  esd->GetRunNumber()));
369
370     Print();
371     fFirstEvent = false;
372   }
373
374   // Some variables 
375   UInt_t   triggers; // Trigger bits
376   Bool_t   lowFlux;  // Low flux flag
377   UShort_t iVz;      // Vertex bin from ESD
378   Double_t vZ;       // Z coordinate from ESD
379   Double_t cent;     // Centrality 
380   UShort_t iVzMc;    // Vertex bin from MC
381   Double_t vZMc;     // Z coordinate of IP vertex from MC
382   Double_t b;        // Impact parameter
383   Int_t    nPart;    // Number of participants 
384   Int_t    nBin;     // Number of binary collisions 
385   Double_t phiR;     // Reaction plane from MC
386   UShort_t nClusters;// Number of SPD clusters 
387   // Process the data 
388   UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, 
389                                      cent, nClusters);
390   fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, 
391                        b, nPart, nBin, phiR);
392
393   Bool_t isInel   = triggers & AliAODForwardMult::kInel;
394   Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
395
396   // Fill the event count histograms 
397   if (isInel)           fHEventsTr->Fill(vZMc);
398   if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
399   fHEvents->Fill(vZMc);
400
401   // Now find our vertex bin object 
402   VtxBin* bin = 0;
403   if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins()) 
404     bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
405   if (!bin) { 
406     // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
407     return;
408   }
409
410   // Clear our ESD object 
411   fESDFMD.Clear();
412
413   // Get FMD data 
414   AliESDFMD* esdFMD = esd->GetFMDData();
415   
416   // Now process our input data and store in own ESD object 
417   fTrackDensity.Calculate(*esdFMD, *mcEvent, vZMc, fESDFMD, bin->fPrimary);
418   bin->fCounts->Fill(0.5);
419
420   // And then bin the data in our vtxbin 
421   for (UShort_t d=1; d<=3; d++) { 
422     UShort_t nr = (d == 1 ? 1 : 2);
423     for (UShort_t q=0; q<nr; q++) { 
424       Char_t      r = (q == 0 ? 'I' : 'O');
425       UShort_t    ns= (q == 0 ?  20 :  40);
426       UShort_t    nt= (q == 0 ? 512 : 256);
427       TH2D*       h = bin->fHists.Get(d,r);
428
429       for (UShort_t s=0; s<ns; s++) { 
430         for (UShort_t t=0; t<nt; t++) {
431           Float_t mult = fESDFMD.Multiplicity(d,r,s,t);
432           
433           if (mult == 0 || mult > 20) continue;
434
435           Float_t phi = fESDFMD.Phi(d,r,s,t) / 180 * TMath::Pi();
436           Float_t eta = fESDFMD.Eta(d,r,s,t);
437           h->Fill(eta,phi,mult);
438         } // for t
439       } // for s 
440     } // for q 
441   } // for d
442 }
443
444 //____________________________________________________________________
445 void
446 AliForwardMCCorrectionsTask::Terminate(Option_t*)
447 {
448   // 
449   // End of job
450   // 
451   // Parameters:
452   //    option Not used 
453   //
454
455   fList = dynamic_cast<TList*>(GetOutputData(1));
456   if (!fList) {
457     AliError("No output list defined");
458     return;
459   }
460
461   // Output list 
462   TList* output = new TList;
463   output->SetOwner();
464   output->SetName(Form("%sResults", GetName()));
465
466   // --- Fill correction object --------------------------------------
467   AliFMDCorrSecondaryMap* corr = new AliFMDCorrSecondaryMap;
468   corr->SetVertexAxis(fVtxAxis);
469   corr->SetEtaAxis(fEtaAxis);
470   
471   TIter     next(fVtxBins);
472   VtxBin*   bin = 0;
473   UShort_t  iVz = 1;
474   while ((bin = static_cast<VtxBin*>(next()))) 
475     bin->Finish(fList, output, iVz++, corr);
476
477   output->Add(corr);
478
479   PostData(2, output);
480 }
481
482 //____________________________________________________________________
483 void
484 AliForwardMCCorrectionsTask::Print(Option_t* option) const
485 {
486   std::cout << ClassName() << "\n"
487             << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
488             << "  Vertex range:     [" << fVtxAxis.GetXmin() 
489             << "," << fVtxAxis.GetXmax() << "]\n"
490             << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
491             << "  Eta range:        [" << fEtaAxis.GetXmin() 
492             << "," << fEtaAxis.GetXmax() << "]"
493             << std::endl;
494   gROOT->IncreaseDirLevel();
495   fInspector.Print(option);
496   fTrackDensity.Print(option);
497   gROOT->DecreaseDirLevel();
498 }
499
500 //====================================================================
501 const char*
502 AliForwardMCCorrectionsTask::VtxBin::BinName(Double_t low, 
503                                              Double_t high) 
504 {
505   static TString buf;
506   buf = Form("vtx%+05.1f_%+05.1f", low, high);
507   buf.ReplaceAll("+", "p");
508   buf.ReplaceAll("-", "m");
509   buf.ReplaceAll(".", "d");
510   return buf.Data();
511 }
512
513
514 //____________________________________________________________________
515 AliForwardMCCorrectionsTask::VtxBin::VtxBin()
516   : fHists(), 
517     fPrimary(0),
518     fCounts(0)
519 {
520 }
521 //____________________________________________________________________
522 AliForwardMCCorrectionsTask::VtxBin::VtxBin(Double_t low, 
523                                             Double_t high, 
524                                             const TAxis& axis)
525   : TNamed(BinName(low, high), 
526            Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
527     fHists(), 
528     fPrimary(0),
529     fCounts(0)
530 {
531   fHists.Init(axis);
532
533   fPrimary = new TH2D("primary", "Primaries", 
534                       axis.GetNbins(), axis.GetXmin(), axis.GetXmax(), 
535                       40, 0, 2*TMath::Pi());
536   fPrimary->SetXTitle("#eta");
537   fPrimary->SetYTitle("#varphi [radians]");
538   fPrimary->Sumw2();
539   fPrimary->SetDirectory(0);
540
541   fCounts = new TH1D("counts", "Counts", 1, 0, 1);
542   fCounts->SetXTitle("Events");
543   fCounts->SetYTitle("# of Events");
544   fCounts->SetDirectory(0);
545 }
546
547 //____________________________________________________________________
548 AliForwardMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
549   : TNamed(o),
550     fHists(o.fHists),
551     fPrimary(0), 
552     fCounts(0)
553 {
554   if (o.fPrimary) {
555     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
556     fPrimary->SetDirectory(0);
557   }
558   if (o.fCounts) {
559     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
560     fCounts->SetDirectory(0);
561   }
562 }
563
564 //____________________________________________________________________
565 AliForwardMCCorrectionsTask::VtxBin&
566 AliForwardMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
567 {
568   TNamed::operator=(o);
569   fHists   = o.fHists;
570   fPrimary = 0;
571   fCounts  = 0;
572   if (o.fPrimary) {
573     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
574     fPrimary->SetDirectory(0);
575   }
576   if (o.fCounts) {
577     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
578     fCounts->SetDirectory(0);
579   }
580   return *this;
581 }
582
583 //____________________________________________________________________
584 void
585 AliForwardMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
586 {
587   TList* d = new TList;
588   d->SetName(GetName());
589   d->SetOwner();
590   l->Add(d);
591
592   d->Add(fHists.fFMD1i);
593   d->Add(fHists.fFMD2i);
594   d->Add(fHists.fFMD2o);
595   d->Add(fHists.fFMD3i);
596   d->Add(fHists.fFMD3o);
597
598   d->Add(fPrimary);
599   d->Add(fCounts);
600 }
601
602 //____________________________________________________________________
603 TH2D*
604 AliForwardMCCorrectionsTask::VtxBin::MakeBg(const TH2D* hits, 
605                                             const TH2D* primary) const
606 {
607   TH2D* h = static_cast<TH2D*>(hits->Clone());
608   h->SetDirectory(0);
609   TString n(h->GetName());
610   n.ReplaceAll("_cache", "");
611   h->SetName(n);
612   h->Divide(primary);
613   
614   return h;
615 }
616   
617 //____________________________________________________________________
618 void
619 AliForwardMCCorrectionsTask::VtxBin::Finish(const TList* input, 
620                                             TList* output, 
621                                             UShort_t iVz,
622                                             AliFMDCorrSecondaryMap* map)
623 {
624   TList* out = new TList;
625   out->SetName(GetName());
626   out->SetOwner();
627   output->Add(out);
628
629   TList* l = static_cast<TList*>(input->FindObject(GetName()));
630   if (!l) { 
631     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
632     return;
633   }
634
635   TH2D*   fmd1i = static_cast<TH2D*>(l->FindObject("FMD1I_cache"));
636   TH2D*   fmd2i = static_cast<TH2D*>(l->FindObject("FMD2I_cache"));
637   TH2D*   fmd2o = static_cast<TH2D*>(l->FindObject("FMD2O_cache"));
638   TH2D*   fmd3i = static_cast<TH2D*>(l->FindObject("FMD3I_cache"));
639   TH2D*   fmd3o = static_cast<TH2D*>(l->FindObject("FMD3O_cache"));
640   TH2D*   primO = static_cast<TH2D*>(l->FindObject("primary"));
641   if (!fmd1i || !fmd2i || !fmd2o || !fmd3i || !fmd3o || !primO) {
642     AliError(Form("Missing histogram(s): %p,%p,%p,%p,%p,%p",
643                   fmd1i, fmd2i, fmd2o, fmd3i, fmd3o, primO));
644     return;
645   }
646
647   // Half coverage in phi for inners
648   TH2D*   primI = static_cast<TH2D*>(primO->Clone());
649   primI->SetDirectory(0);
650   primI->RebinY(2); 
651
652   TH2D* bg1i = MakeBg(fmd1i, primI);
653   TH2D* bg2i = MakeBg(fmd2i, primI);
654   TH2D* bg2o = MakeBg(fmd2o, primO);
655   TH2D* bg3i = MakeBg(fmd3i, primI);
656   TH2D* bg3o = MakeBg(fmd3o, primO);
657   map->SetCorrection(1, 'I', iVz, bg1i);
658   map->SetCorrection(2, 'I', iVz, bg2i);
659   map->SetCorrection(2, 'O', iVz, bg2o);
660   map->SetCorrection(3, 'I', iVz, bg3i);
661   map->SetCorrection(3, 'O', iVz, bg3o);
662   out->Add(bg1i);
663   out->Add(bg2i);
664   out->Add(bg2o);
665   out->Add(bg3i);
666   out->Add(bg3o);
667  
668 }
669
670 //
671 // EOF
672 //