Various improvements
[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 "AliFMDStripIndex.h"
27 #include <TH1.h>
28 #include <TH3D.h>
29 #include <TDirectory.h>
30 #include <TTree.h>
31
32 //====================================================================
33 namespace {
34   const char* GetEventName(Bool_t tr, Bool_t vtx) 
35   {
36     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
37   }
38   const char* GetHitsName(UShort_t d, Char_t r) 
39   {
40     return Form("hitsFMD%d%c", d, r);
41   }
42   const char* GetStripsName(UShort_t d, Char_t r)
43   {
44     return Form("stripsFMD%d%c", d, r);
45   }
46   const char* GetPrimaryName(Char_t r, Bool_t trVtx)
47   {
48     return Form("primaries%s%s", 
49                 ((r == 'I' || r == 'i') ? "Inner" : "Outer"), 
50                 (trVtx ? "TrVtx" : "All"));
51   }
52 }
53
54 //====================================================================
55 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
56   : AliAnalysisTaskSE(),
57     fHEvents(0), 
58     fHEventsTr(0), 
59     fHEventsTrVtx(0),
60     fHEventsVtx(0), 
61     fHTriggers(0),
62     fPrimaryInnerAll(0),   
63     fPrimaryOuterAll(0),   
64     fPrimaryInnerTrVtx(0), 
65     fPrimaryOuterTrVtx(0), 
66     fHitsFMD1i(0),         
67     fHitsFMD2i(0),         
68     fHitsFMD2o(0),         
69     fHitsFMD3i(0),         
70     fHitsFMD3o(0),         
71     fStripsFMD1i(0),       
72     fStripsFMD2i(0),       
73     fStripsFMD2o(0),       
74     fStripsFMD3i(0),       
75     fStripsFMD3o(0),       
76     fVtxAxis(),
77     fEtaAxis(),
78     fList()
79 {
80   // 
81   // Constructor 
82   // 
83   // Parameters:
84   //    name Name of task 
85   //
86 }
87
88 //____________________________________________________________________
89 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
90   : AliAnalysisTaskSE(name), 
91     fHEvents(0), 
92     fHEventsTr(0), 
93     fHEventsTrVtx(0),
94     fHEventsVtx(0), 
95     fHTriggers(0),
96     fPrimaryInnerAll(0),   
97     fPrimaryOuterAll(0),   
98     fPrimaryInnerTrVtx(0), 
99     fPrimaryOuterTrVtx(0), 
100     fHitsFMD1i(0),         
101     fHitsFMD2i(0),         
102     fHitsFMD2o(0),         
103     fHitsFMD3i(0),         
104     fHitsFMD3o(0),         
105     fStripsFMD1i(0),       
106     fStripsFMD2i(0),       
107     fStripsFMD2o(0),       
108     fStripsFMD3i(0),       
109     fStripsFMD3o(0),       
110     fVtxAxis(10,-10,10), 
111     fEtaAxis(200,-4,6),
112     fList()
113 {
114   // 
115   // Constructor 
116   // 
117   // Parameters:
118   //    name Name of task 
119   //
120   DefineOutput(1, TList::Class());
121 }
122
123 //____________________________________________________________________
124 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
125   : AliAnalysisTaskSE(o),
126     fHEvents(o.fHEvents), 
127     fHEventsTr(o.fHEventsTr), 
128     fHEventsTrVtx(o.fHEventsTrVtx),
129     fHEventsVtx(o.fHEventsVtx), 
130     fHTriggers(o.fHTriggers),
131     fPrimaryInnerAll(o.fPrimaryInnerAll),   
132     fPrimaryOuterAll(o.fPrimaryOuterAll),   
133     fPrimaryInnerTrVtx(o.fPrimaryInnerTrVtx), 
134     fPrimaryOuterTrVtx(o.fPrimaryOuterTrVtx), 
135     fHitsFMD1i(o.fHitsFMD1i),         
136     fHitsFMD2i(o.fHitsFMD2i),         
137     fHitsFMD2o(o.fHitsFMD2o),         
138     fHitsFMD3i(o.fHitsFMD3i),         
139     fHitsFMD3o(o.fHitsFMD3o),         
140     fStripsFMD1i(o.fStripsFMD1i),       
141     fStripsFMD2i(o.fStripsFMD2i),       
142     fStripsFMD2o(o.fStripsFMD2o),       
143     fStripsFMD3i(o.fStripsFMD3i),       
144     fStripsFMD3o(o.fStripsFMD3o),       
145     fVtxAxis(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), o.fVtxAxis.GetXmax()),
146     fEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax()),
147     fList(o.fList)
148 {
149   // 
150   // Copy constructor 
151   // 
152   // Parameters:
153   //    o Object to copy from 
154   //
155 }
156
157 //____________________________________________________________________
158 AliForwardMCCorrectionsTask&
159 AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
160 {
161   // 
162   // Assignment operator 
163   // 
164   // Parameters:
165   //    o Object to assign from 
166   // 
167   // Return:
168   //    Reference to this object 
169   //
170   fHEventsTr         = o.fHEventsTr;
171   fHEventsTrVtx      = o.fHEventsTrVtx;
172   fHTriggers         = o.fHTriggers;
173   SetVertexAxis(o.fVtxAxis);
174   SetEtaAxis(o.fEtaAxis);
175
176   return *this;
177 }
178
179 //____________________________________________________________________
180 void
181 AliForwardMCCorrectionsTask::Init()
182 {
183   // 
184   // Initialize the task 
185   // 
186   //
187 }
188
189 //____________________________________________________________________
190 void
191 AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min, 
192                                            Double_t max)
193 {
194   // 
195   // Set the vertex axis to use
196   // 
197   // Parameters:
198   //    nBins Number of bins
199   //    vzMin Least @f$z@f$ coordinate of interation point
200   //    vzMax Largest @f$z@f$ coordinate of interation point
201   //
202   if (max < min) max = -min;
203   if (min < max) { 
204     Double_t tmp = min;
205     min          = max;
206     max          = tmp;
207   }
208   if (min < -10) 
209     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
210   if (max > +10) 
211     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
212   fVtxAxis.Set(nBin, min, max);
213 }
214 //____________________________________________________________________
215 void
216 AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
217 {
218   // 
219   // Set the vertex axis to use
220   // 
221   // Parameters:
222   //    axis Axis
223   //
224   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
225 }
226
227 //____________________________________________________________________
228 void
229 AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
230 {
231   // 
232   // Set the eta axis to use
233   // 
234   // Parameters:
235   //    nBins Number of bins
236   //    vzMin Least @f$\eta@f$ 
237   //    vzMax Largest @f$\eta@f$ 
238   //
239   if (max < min) max = -min;
240   if (min < max) { 
241     Double_t tmp = min;
242     min          = max;
243     max          = tmp;
244   }
245   if (min < -4) 
246     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
247   if (max > +6) 
248     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
249   fVtxAxis.Set(nBin, min, max);
250 }
251 //____________________________________________________________________
252 void
253 AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
254 {
255   // 
256   // Set the eta axis to use
257   // 
258   // Parameters:
259   //    axis Axis
260   //
261   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
262 }
263   
264 //____________________________________________________________________
265 TH3D*
266 AliForwardMCCorrectionsTask::Make3D(const char* name, const char* title,
267                                Int_t nPhi) const
268 {
269   // 
270   // Make a 3D histogram
271   // 
272   // Parameters:
273   //    name   Name 
274   //    title  Title 
275   //    nPhi   Number of phi bins
276   // 
277   // Return:
278   //    Histogram
279   //
280   TH3D* ret = new TH3D(name, title,
281                        fVtxAxis.GetNbins(), 
282                        fVtxAxis.GetXmin(), 
283                        fVtxAxis.GetXmax(), 
284                        fEtaAxis.GetNbins(), 
285                        fEtaAxis.GetXmin(), 
286                        fEtaAxis.GetXmax(), 
287                        nPhi, 0, 2*TMath::Pi());
288   ret->SetXTitle("v_{z} [cm]");
289   ret->SetYTitle("#eta");
290   ret->SetZTitle("#varphi [radians]");
291   ret->SetDirectory(0);
292   ret->SetStats(0);
293   ret->Sumw2();
294
295   return ret;
296 }
297 //____________________________________________________________________
298 TH1D*
299 AliForwardMCCorrectionsTask::Make1D(const char* name, const char* title) const
300 {
301   // 
302   // Make 1D histogram
303   // 
304   // Parameters:
305   //    name   Name 
306   //    title  Title
307   // 
308   // Return:
309   //    Histogram
310   //
311   TH1D* ret = new TH1D(name, title,
312                        fEtaAxis.GetNbins(), 
313                        fEtaAxis.GetXmin(), 
314                        fEtaAxis.GetXmax());
315   ret->SetXTitle("#eta");
316   ret->SetFillColor(kRed+1);
317   ret->SetFillStyle(3001);
318   ret->SetDirectory(0);
319   ret->SetStats(0);
320   ret->Sumw2();
321
322   return ret;
323 }
324
325 //____________________________________________________________________
326 void
327 AliForwardMCCorrectionsTask::UserCreateOutputObjects()
328 {
329   // 
330   // Create output objects 
331   // 
332   //
333   fList = new TList;
334   fList->SetName(GetName());
335
336   fHEvents = new TH1I(GetEventName(false,false),
337                       "Number of all events", 
338                       fVtxAxis.GetNbins(), 
339                       fVtxAxis.GetXmin(), 
340                       fVtxAxis.GetXmax());
341   fHEvents->SetXTitle("v_{z} [cm]");
342   fHEvents->SetYTitle("# of events");
343   fHEvents->SetFillColor(kBlue+1);
344   fHEvents->SetFillStyle(3001);
345   fHEvents->SetDirectory(0);
346   fList->Add(fHEvents);
347
348   fHEventsTr = new TH1I(GetEventName(true, false), 
349                         "Number of triggered events",
350                         fVtxAxis.GetNbins(), 
351                         fVtxAxis.GetXmin(), 
352                         fVtxAxis.GetXmax());
353   fHEventsTr->SetXTitle("v_{z} [cm]");
354   fHEventsTr->SetYTitle("# of events");
355   fHEventsTr->SetFillColor(kRed+1);
356   fHEventsTr->SetFillStyle(3001);
357   fHEventsTr->SetDirectory(0);
358   fList->Add(fHEventsTr);
359
360   fHEventsTrVtx = new TH1I(GetEventName(true,true),
361                            "Number of events w/trigger and vertex", 
362                            fVtxAxis.GetNbins(), 
363                            fVtxAxis.GetXmin(), 
364                            fVtxAxis.GetXmax());
365   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
366   fHEventsTrVtx->SetYTitle("# of events");
367   fHEventsTrVtx->SetFillColor(kBlue+1);
368   fHEventsTrVtx->SetFillStyle(3001);
369   fHEventsTrVtx->SetDirectory(0);
370   fList->Add(fHEventsTrVtx);
371   
372   fHEventsVtx = new TH1I(GetEventName(false,true),
373                          "Number of events w/vertex", 
374                          fVtxAxis.GetNbins(), 
375                          fVtxAxis.GetXmin(), 
376                          fVtxAxis.GetXmax());
377   fHEventsVtx->SetXTitle("v_{z} [cm]");
378   fHEventsVtx->SetYTitle("# of events");
379   fHEventsVtx->SetFillColor(kBlue+1);
380   fHEventsVtx->SetFillStyle(3001);
381   fHEventsVtx->SetDirectory(0);
382   fList->Add(fHEventsVtx);
383
384       
385   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
386   fHTriggers->SetFillColor(kRed+1);
387   fHTriggers->SetFillStyle(3001);
388   fHTriggers->SetStats(0);
389   fHTriggers->SetDirectory(0);
390   fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
391   fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
392   fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
393   fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
394   fHTriggers->GetXaxis()->SetBinLabel(5,"A");
395   fHTriggers->GetXaxis()->SetBinLabel(6,"B");
396   fHTriggers->GetXaxis()->SetBinLabel(7,"C");
397   fHTriggers->GetXaxis()->SetBinLabel(8,"E");
398   fList->Add(fHTriggers);
399
400   fPrimaryInnerAll   = Make3D(GetPrimaryName('I',false), "Primary particles, "
401                               "20 #varphi bins, all events", 20);
402   fPrimaryOuterAll   = Make3D(GetPrimaryName('O',false), "Primary particles, "
403                               "40 #varphi bins, all events", 40);
404   fPrimaryInnerTrVtx = Make3D(GetPrimaryName('I',true), "Primary particles, "
405                               "20 #varphi bins, Tr+Vtx events", 20);
406   fPrimaryOuterTrVtx = Make3D(GetPrimaryName('O',true), "Primary particles, "
407                               "40 #varphi bins, Tr+Vtx events", 40);
408   TList* primaries = new TList;
409   primaries->SetName("primaries");
410   primaries->Add(fPrimaryInnerAll);
411   primaries->Add(fPrimaryOuterAll);
412   primaries->Add(fPrimaryInnerTrVtx);
413   primaries->Add(fPrimaryOuterTrVtx);
414   fList->Add(primaries);
415
416
417   fHitsFMD1i   = Make3D(GetHitsName(1,'i'),   "All hits in FMD1i, Tr+Vtx", 20);
418   fHitsFMD2i   = Make3D(GetHitsName(2,'i'),   "All hits in FMD2i, Tr+Vtx", 20);
419   fHitsFMD2o   = Make3D(GetHitsName(2,'o'),   "All hits in FMD2o, Tr+Vtx", 40);
420   fHitsFMD3i   = Make3D(GetHitsName(3,'i'),   "All hits in FMD3i, Tr+Vtx", 20);
421   fHitsFMD3o   = Make3D(GetHitsName(3,'o'),   "All hits in FMD3o, Tr+Vtx", 40);
422   TList* hits = new TList;
423   hits->SetName("hits");
424   hits->Add(fHitsFMD1i);
425   hits->Add(fHitsFMD2i);
426   hits->Add(fHitsFMD2o);
427   hits->Add(fHitsFMD3i);
428   hits->Add(fHitsFMD3o);
429   fList->Add(hits);
430
431   fStripsFMD1i = Make1D(GetStripsName(1,'i'), "# hit strips in FMD1i, Tr+Vtx");
432   fStripsFMD2i = Make1D(GetStripsName(2,'i'), "# hit strips in FMD2i, Tr+Vtx");
433   fStripsFMD2o = Make1D(GetStripsName(2,'o'), "# hit strips in FMD2o, Tr+Vtx");
434   fStripsFMD3i = Make1D(GetStripsName(3,'i'), "# hit strips in FMD3i, Tr+Vtx");
435   fStripsFMD3o = Make1D(GetStripsName(3,'o'), "# hit strips in FMD3o, Tr+Vtx");
436   TList* strips = new TList;
437   strips->SetName("strips");
438   strips->Add(fStripsFMD1i);
439   strips->Add(fStripsFMD2i);
440   strips->Add(fStripsFMD2o);
441   strips->Add(fStripsFMD3i);
442   strips->Add(fStripsFMD3o);
443   fList->Add(strips);
444
445   PostData(1, fList);
446 }
447 //____________________________________________________________________
448 void
449 AliForwardMCCorrectionsTask::UserExec(Option_t*)
450 {
451   // 
452   // Process each event 
453   // 
454   // Parameters:
455   //    option Not used
456   //  
457
458   // Get the input data - MC event
459   AliMCEvent*  mcEvent = MCEvent();
460   if (!mcEvent) { 
461     AliWarning("No MC event found");
462     return;
463   }
464
465   // Get the input data - ESD event
466   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
467   if (!esd) { 
468     AliWarning("No ESD event found for input event");
469     return;
470   }
471   
472   // Get the particle stack 
473   AliStack* stack = mcEvent->Stack();
474
475   // Get the event generate header 
476   AliHeader*          mcHeader  = mcEvent->Header();
477   AliGenEventHeader*  genHeader = mcHeader->GenEventHeader();
478   
479   // Get the generator vertex 
480   TArrayF mcVertex;
481   genHeader->PrimaryVertex(mcVertex);
482   Double_t mcVtxZ = mcVertex.At(2);
483
484   // Check the MC vertex is in range 
485   Int_t mcVtxBin = fVtxAxis.FindBin(mcVtxZ);
486   if (mcVtxBin < 1 || mcVtxBin > fVtxAxis.GetNbins()) {
487 #ifdef VERBOSE 
488     AliWarning(Form("MC v_z=%f is out of rante [%,%f]", 
489                     mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
490 #endif
491     return;
492   }
493
494   UInt_t   triggers;
495   Bool_t   gotTrigggers;
496   Bool_t   gotInel;
497   Double_t vZ;
498   Bool_t   gotVertex;
499 #if 0
500   // Use event inspector instead 
501   // Get the triggers 
502   UInt_t triggers = 0;
503   Bool_t gotTriggers = AliForwardUtil::ReadTriggers(esd,triggers,fHTriggers);
504   Bool_t gotInel     = triggers & AliAODForwardMult::kInel;
505   
506   // Get the ESD vertex 
507   Double_t vZ = -1000000;
508   Bool_t gotVertex = AliForwardUtil::ReadVertex(esd,vZ);
509 #endif
510
511
512   // Fill the event count histograms 
513   if (gotInel)              fHEventsTr->Fill(mcVtxZ);
514   if (gotInel && gotVertex) fHEventsTrVtx->Fill(mcVtxZ);
515   if (gotVertex)            fHEventsVtx->Fill(mcVtxZ);
516   fHEvents->Fill(mcVtxZ);
517
518   // Cache of the hits 
519   AliFMDFloatMap hitsByStrip;
520   AliFMDFloatMap lastByStrip;
521
522   // Loop over all tracks 
523   Int_t nTracks = mcEvent->GetNumberOfTracks();
524   for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
525     AliMCParticle* particle = 
526       static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
527     
528     // Check the returned particle 
529     if (!particle) continue;
530
531     // Check if this charged and a primary 
532     Bool_t isCharged = particle->Charge() != 0;
533     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
534
535     // Fill (eta,phi) of the particle into histograsm for b
536     Double_t eta = particle->Eta();
537     Double_t phi = particle->Phi();
538     
539     if (isCharged && isPrimary) 
540       FillPrimary(gotInel, gotVertex, mcVtxZ, eta, phi);
541     
542     // For the rest, ignore non-collisions, and events out of vtx range 
543     if (!gotInel || gotVertex) continue;
544     
545     Int_t nTrRef = particle->GetNumberOfTrackReferences();
546     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
547       AliTrackReference* trackRef = particle->GetTrackReference(iTrRef);
548       
549       // Check existence 
550       if (!trackRef) continue;
551
552       // Check that we hit an FMD element 
553       if (trackRef->DetectorId() != AliTrackReference::kFMD) 
554         continue;
555
556       // Get the detector coordinates 
557       UShort_t d, s, t;
558       Char_t r;
559       AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t);
560       
561       // Check if mother (?) is charged and that we haven't dealt with it 
562       // already
563       Int_t lastTrack = Int_t(lastByStrip(d,r,s,t));
564       if (!isCharged || iTr == lastTrack) continue;
565
566       // Increase counter of hits per strip 
567       hitsByStrip(d,r,s,t) += 1;
568
569       // Double_t trRefEta = esd->GetFMDData()->Eta(d,r,s,t);
570       // Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
571
572       // Fill strip information into histograms 
573       FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
574
575       // Set the last processed track number - marking it as done for
576       // this strip
577       lastByStrip(d,r,s,t) = Float_t(iTr);
578       
579       // Flag neighboring strips too 
580       Int_t nStrip = (r == 'I' || r == 'i' ? 512 : 256);
581       if (t > 0)        lastByStrip(d,r,s,t-1) = Float_t(iTr);
582       if (t < nStrip-1) lastByStrip(d,r,s,t+1) = Float_t(iTr);
583     }
584   }
585 }
586
587 //____________________________________________________________________
588 void
589 AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx, 
590                                     Double_t vz, Double_t eta, Double_t phi) 
591 {
592   // 
593   // Fill in primary information
594   // 
595   // Parameters:
596   //    gotInel   Got INEL trigger from ESD
597   //    gotVtx    Got vertex Z from ESD 
598   //    vz        @f$z@f$ coordinate of interation point
599   //    eta       Pseudo rapidity 
600   //    phi       Azimuthal angle
601   //
602   if (gotInel && gotVtx) {
603     // This takes the place of hPrimary_FMD_<r>_vtx<v> and 
604     // Analysed_FMD<r>_vtx<v>
605     fPrimaryInnerTrVtx->Fill(vz,eta,phi);
606     fPrimaryOuterTrVtx->Fill(vz,eta,phi);
607   }
608   // This takes the place of Inel_FMD<r>_vtx<v> and 
609   // Analysed_FMD<r>_vtx<v>
610   fPrimaryInnerAll->Fill(vz,eta,phi);
611   fPrimaryOuterAll->Fill(vz,eta,phi);
612 }
613
614 //____________________________________________________________________
615 void
616 AliForwardMCCorrectionsTask::FillStrip(UShort_t d, Char_t r, 
617                                   Double_t vz, Double_t eta, Double_t phi,
618                                   Bool_t first) 
619 {
620   // 
621   // Fill in per-strip information
622   // 
623   // Parameters:
624   //    d         Detector
625   //    r         Ring
626   //    vz        @f$z@f$ coordinate of interation point
627   //    eta       Pseudo rapidity 
628   //    phi       Azimuthal angle
629   //    first     First fill in this event
630   //
631
632   // Number of hit strips per eta bin 
633   TH1D* strips = 0; 
634   // All hits per ring, vertex in eta,phi bins.  This takes the place
635   // of hHits_FMD<d><r>_vtx<v> and DoubleHits_FMD<d><r> (the later
636   // being just the 2D projection over the X bins)
637   TH3D* hits   = 0; 
638   switch (d) { 
639   case 1: 
640     hits   = fHitsFMD1i; 
641     strips = fStripsFMD1i; 
642     break;
643   case 2: 
644     hits   = (r == 'I' || r == 'i' ? fHitsFMD2i   : fHitsFMD2o);
645     strips = (r == 'I' || r == 'i' ? fStripsFMD2i : fStripsFMD2o);
646     break;
647   case 3: 
648     hits   = (r == 'I' || r == 'i' ? fHitsFMD3i   : fHitsFMD3o);
649     strips = (r == 'I' || r == 'i' ? fStripsFMD3i : fStripsFMD3o);
650     break;
651   }
652   if (!hits || !strips) return;
653   
654   if (first) strips->Fill(eta);
655   
656   hits->Fill(vz, eta, phi);
657 }
658
659 //____________________________________________________________________
660 TH2D*
661 AliForwardMCCorrectionsTask::GetVertexProj(Int_t v, TH3D* src) const
662 {
663   // 
664   // Get vertex project
665   // 
666   // Parameters:
667   //    v   Vertex bin 
668   //    src Source 3D histogram 
669   // 
670   // Return:
671   //    2D projection of the V'th bin
672   //
673   Int_t nX = src->GetNbinsX();
674   if (v > nX || v < 1) return 0;
675
676   src->GetXaxis()->SetRange(v,v+1);
677   
678   TH2D* ret   = static_cast<TH2D*>(src->Project3D("zy"));
679   ret->SetName(Form("%s_vtx%02d", src->GetName(), v));
680   ret->SetDirectory(0);
681
682   src->GetXaxis()->SetRange(0,nX+1);
683
684   return ret;
685 }
686
687
688 //____________________________________________________________________
689 void
690 AliForwardMCCorrectionsTask::Terminate(Option_t*)
691 {
692   // 
693   // End of job
694   // 
695   // Parameters:
696   //    option Not used 
697   //
698
699   TList* list = dynamic_cast<TList*>(GetOutputData(1));
700   if (!list) {
701     AliError("No output list defined");
702     return;
703   }
704
705   TList* primaries = static_cast<TList*>(list->FindObject("primaries"));
706   TList* hits      = static_cast<TList*>(list->FindObject("hits"));
707   TList* strips    = static_cast<TList*>(list->FindObject("strips"));
708   if (!primaries || !hits || !strips) { 
709     AliError(Form("Could not find all sub-lists in %s (p=%p,h=%p,s=%p)",
710                   list->GetName(), primaries, hits, strips));
711     return;
712   }
713
714   TH1I* eventsAll = 
715     static_cast<TH1I*>(list->FindObject(GetEventName(false,false)));
716   TH1I* eventsTr = 
717     static_cast<TH1I*>(list->FindObject(GetEventName(true,false)));
718   TH1I* eventsVtx = 
719     static_cast<TH1I*>(list->FindObject(GetEventName(false,true)));
720   TH1I* eventsTrVtx = 
721     static_cast<TH1I*>(list->FindObject(GetEventName(true,true)));
722   if (!eventsAll || !eventsTr || !eventsVtx || !eventsTrVtx) {
723     AliError(Form("Could not find all event histograms in %s "
724                   "(a=%p,t=%p,v=%p,tv=%p)", list->GetName(), 
725                   eventsAll, eventsTr, eventsVtx, eventsTrVtx));
726     return;
727   }
728
729   TH3D* primIAll = 
730     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('I',false)));
731   TH3D* primOAll = 
732     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',false)));
733   TH3D* primITrVtx = 
734     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('I',true)));
735   TH3D* primOTrVtx = 
736     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',true)));
737   if (!primIAll || !primOAll || !primITrVtx || !primOTrVtx) {
738     AliError(Form("Could not find all primary particle histograms in %s "
739                   "(ai=%p,ao=%p,tvi=%p,tvo=%p)", list->GetName(), 
740                   primIAll, primOAll, primITrVtx, primOTrVtx));
741     return;
742   }
743     
744   TH3D* hits1i = static_cast<TH3D*>(hits->FindObject(GetHitsName(1,'i')));
745   TH3D* hits2i = static_cast<TH3D*>(hits->FindObject(GetHitsName(2,'i')));
746   TH3D* hits2o = static_cast<TH3D*>(hits->FindObject(GetHitsName(2,'o')));
747   TH3D* hits3i = static_cast<TH3D*>(hits->FindObject(GetHitsName(3,'i')));
748   TH3D* hits3o = static_cast<TH3D*>(hits->FindObject(GetHitsName(3,'o')));
749   if (!hits1i || !hits2i || !hits2o || hits3i || hits3o) {
750     AliError(Form("Could not find all ring hits histograms in %s " 
751                   "(1i=%p,2i=%p,2o=%p,3i=%p,3o=%p)", hits->GetName(), 
752                   hits1i, hits2i, hits2o, hits3i, hits3o));
753     return;
754   }
755
756   TH1D* strips1i = static_cast<TH1D*>(strips->FindObject(GetStripsName(1,'i')));
757   TH1D* strips2i = static_cast<TH1D*>(strips->FindObject(GetStripsName(2,'i')));
758   TH1D* strips2o = static_cast<TH1D*>(strips->FindObject(GetStripsName(2,'o')));
759   TH1D* strips3i = static_cast<TH1D*>(strips->FindObject(GetStripsName(3,'i')));
760   TH1D* strips3o = static_cast<TH1D*>(strips->FindObject(GetStripsName(3,'o')));
761   if (!strips1i || !strips2i || !strips2o || strips3i || strips3o) {
762     AliError(Form("Could not find all ring strips histograms in %s " 
763                   "(1i=%p,2i=%p,2o=%p,3i=%p,3o=%p)", strips->GetName(), 
764                   strips1i, strips2i, strips2o, strips3i, strips3o));
765     return;
766   }
767
768   // Calculate the over-all event selection efficiency 
769   TH1D* selEff = new TH1D("selEff", "Event selection efficiency", 
770                           fVtxAxis.GetNbins(), 
771                           fVtxAxis.GetXmin(),  
772                           fVtxAxis.GetXmax());
773   selEff->Sumw2();
774   selEff->SetDirectory(0);
775   selEff->SetFillColor(kMagenta+1);
776   selEff->SetFillStyle(3001);
777   selEff->Add(eventsAll);
778   selEff->Divide(eventsTrVtx);
779   list->Add(selEff);
780
781   // Loop over vertex bins and do vertex dependendt stuff 
782   for (Int_t v = 1; v <= fVtxAxis.GetNbins(); v++) {
783     // Make a sub-list in the output 
784     TList* vl = new TList;
785     vl->SetName(Form("vtx%02d", v));
786     list->Add(vl);
787
788     // Get event counts 
789     Int_t nEventsAll   = eventsAll->GetBinContent(v);
790     Int_t nEventsTr    = eventsTr->GetBinContent(v);
791     Int_t nEventsVtx   = eventsVtx->GetBinContent(v);
792     Int_t nEventsTrVtx = eventsTrVtx->GetBinContent(v);
793
794     // Project event histograms, set names, and store  
795     TH2D* primIAllV   = GetVertexProj(v, primIAll);
796     TH2D* primOAllV   = GetVertexProj(v, primOAll);
797     TH2D* primITrVtxV = GetVertexProj(v, primITrVtx);
798     TH2D* primOTrVtxV = GetVertexProj(v, primOTrVtx);
799     vl->Add(primIAllV);
800     vl->Add(primOAllV);
801     vl->Add(primITrVtxV);
802     vl->Add(primOTrVtxV);
803     
804     primIAllV->Scale(1. / nEventsAll);
805     primOAllV->Scale(1. / nEventsAll);
806     primITrVtxV->Scale(1. / nEventsTr);
807     primOTrVtxV->Scale(1. / nEventsTr);
808
809     // Calculate the vertex bias on the d^2N/detadphi 
810     TH2D* selBiasI = 
811       static_cast<TH2D*>(primITrVtxV->Clone(Form("selBiasI%02d",v)));
812     TH2D* selBiasO = 
813       static_cast<TH2D*>(primOTrVtxV->Clone(Form("selBiasO%02d",v)));
814     selBiasI->SetTitle(Form("Event selection bias in vertex bin %d", v));
815     selBiasO->SetTitle(Form("Event selection bias in vertex bin %d", v));
816     selBiasI->Divide(primIAllV);
817     selBiasO->Divide(primOAllV);
818     vl->Add(selBiasI);
819     vl->Add(selBiasO);
820
821     for (UShort_t d = 1; d <= 3; d++) { 
822       UShort_t nQ = (d == 1 ? 1 : 2);
823       for (UShort_t q = 0; q < nQ; q++) { 
824         Char_t   r  = (q == 0 ? 'I' : 'O');
825         TH3D* hits3D = 0;
826         TH2D* prim2D = (q == 0 ? primITrVtxV : primOTrVtxV);
827         switch (d) { 
828         case 1: hits3D = hits1i; break;
829         case 2: hits3D = (q == 0 ? hits2i : hits2o); break;
830         case 3: hits3D = (q == 0 ? hits3i : hits3o); break;
831         }
832
833         TH2D* sec = GetVertexProj(v, hits3D);
834         sec->SetName(Form("secondaryFMD%d%c_vtx%02d", d, r, v));
835         sec->SetTitle(Form("Secondary correction map for FMD%d%c "
836                            "in vertex bin %d", d, r, v));
837         sec->Divide(prim2D);
838         vl->Add(sec);
839
840         if (v > 1) continue;
841         
842         // Do the double hit correction (only done once per ring in
843         // the vertex loop)
844         TH1D* hStrips = 0;
845         switch (d) { 
846         case 1: hStrips = strips1i; break;
847         case 2: hStrips = (q == 0 ? strips2i : strips2o); break;
848         case 3: hStrips = (q == 0 ? strips3i : strips3o); break;
849         }
850
851         TH2D* hits2D    = GetVertexProj(v, hits3D);
852         TH1D* doubleHit = hits2D->ProjectionX(Form("doubleHitFMD%d%c",d,r));
853         doubleHit->SetTitle(Form("Double hit correction for FMD%d%c",d,r));
854         doubleHit->SetDirectory(0);
855         doubleHit->SetFillColor(kGreen+1);
856         doubleHit->SetFillStyle(3001);
857         doubleHit->Sumw2();
858         doubleHit->Divide(hStrips);
859         list->Add(doubleHit);
860       }
861     }
862   }
863 }
864
865 //____________________________________________________________________
866 void
867 AliForwardMCCorrectionsTask::Print(Option_t*) const
868 {
869 }
870
871 //
872 // EOF
873 //