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