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