More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMCCorrectionsTask.cxx
CommitLineData
0bd4b00f 1#include "AliForwardMCCorrectionsTask.h"
563a673f 2#include "AliTriggerAnalysis.h"
3#include "AliPhysicsSelection.h"
4#include "AliLog.h"
0bd4b00f 5// #include "AliFMDAnaParameters.h"
6#include "AliHeader.h"
7#include "AliGenEventHeader.h"
563a673f 8#include "AliESDEvent.h"
9#include "AliAODHandler.h"
10#include "AliMultiplicity.h"
11#include "AliInputEventHandler.h"
0bd4b00f 12#include "AliStack.h"
13#include "AliMCEvent.h"
14#include "AliFMDStripIndex.h"
563a673f 15#include <TH1.h>
0bd4b00f 16#include <TH3D.h>
563a673f 17#include <TDirectory.h>
18#include <TTree.h>
19
20//====================================================================
21namespace {
22 const char* GetEventName(Bool_t tr, Bool_t vtx)
23 {
0bd4b00f 24 return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
563a673f 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//====================================================================
0bd4b00f 43AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
563a673f 44 : AliAnalysisTaskSE(),
0bd4b00f 45 fHEvents(0),
563a673f 46 fHEventsTr(0),
47 fHEventsTrVtx(0),
0bd4b00f 48 fHEventsVtx(0),
563a673f 49 fHTriggers(0),
0bd4b00f 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),
563a673f 64 fVtxAxis(),
0bd4b00f 65 fEtaAxis(),
66 fList()
563a673f 67{
68}
69
70//____________________________________________________________________
0bd4b00f 71AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
563a673f 72 : AliAnalysisTaskSE(name),
0bd4b00f 73 fHEvents(0),
563a673f 74 fHEventsTr(0),
0bd4b00f 75 fHEventsTrVtx(0),
76 fHEventsVtx(0),
563a673f 77 fHTriggers(0),
0bd4b00f 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),
563a673f 92 fVtxAxis(10,-10,10),
0bd4b00f 93 fEtaAxis(200,-4,6),
94 fList()
563a673f 95{
96 DefineOutput(1, TList::Class());
97}
98
99//____________________________________________________________________
0bd4b00f 100AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
563a673f 101 : AliAnalysisTaskSE(o),
0bd4b00f 102 fHEvents(o.fHEvents),
563a673f 103 fHEventsTr(o.fHEventsTr),
0bd4b00f 104 fHEventsTrVtx(o.fHEventsTrVtx),
105 fHEventsVtx(o.fHEventsVtx),
563a673f 106 fHTriggers(o.fHTriggers),
0bd4b00f 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)
563a673f 124{
125}
126
127//____________________________________________________________________
0bd4b00f 128AliForwardMCCorrectionsTask&
129AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
563a673f 130{
131 fHEventsTr = o.fHEventsTr;
132 fHEventsTrVtx = o.fHEventsTrVtx;
133 fHTriggers = o.fHTriggers;
0bd4b00f 134 SetVertexAxis(o.fVtxAxis);
135 SetEtaAxis(o.fEtaAxis);
563a673f 136
137 return *this;
138}
139
140//____________________________________________________________________
141void
0bd4b00f 142AliForwardMCCorrectionsTask::Init()
563a673f 143{
144}
145
146//____________________________________________________________________
147void
0bd4b00f 148AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
149 Double_t max)
563a673f 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//____________________________________________________________________
164void
0bd4b00f 165AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
563a673f 166{
167 SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
168}
169
170//____________________________________________________________________
171void
0bd4b00f 172AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
563a673f 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//____________________________________________________________________
187void
0bd4b00f 188AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
563a673f 189{
190 SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
191}
192
563a673f 193//____________________________________________________________________
194TH3D*
0bd4b00f 195AliForwardMCCorrectionsTask::Make3D(const char* name, const char* title,
563a673f 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//____________________________________________________________________
216TH1D*
0bd4b00f 217AliForwardMCCorrectionsTask::Make1D(const char* name, const char* title) const
563a673f 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//____________________________________________________________________
234void
0bd4b00f 235AliForwardMCCorrectionsTask::UserCreateOutputObjects()
563a673f 236{
237 fList = new TList;
238 fList->SetName(GetName());
239
0bd4b00f 240 fHEvents = new TH1I(GetEventName(false,false),
563a673f 241 "Number of all events",
0bd4b00f 242 fVtxAxis.GetNbins(),
243 fVtxAxis.GetXmin(),
244 fVtxAxis.GetXmax());
563a673f 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
0bd4b00f 252 fHEventsTr = new TH1I(GetEventName(true, false),
253 "Number of triggered events",
563a673f 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
0bd4b00f 264 fHEventsTrVtx = new TH1I(GetEventName(true,true),
563a673f 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);
0bd4b00f 275
276 fHEventsVtx = new TH1I(GetEventName(false,true),
563a673f 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//____________________________________________________________________
352void
0bd4b00f 353AliForwardMCCorrectionsTask::UserExec(Option_t*)
563a673f 354{
355 // Get the input data - MC event
356 AliMCEvent* mcEvent = MCEvent();
0bd4b00f 357 if (mcEvent) {
563a673f 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();
0bd4b00f 374 AliGenEventHeader* genHeader = mcHeader->GenEventHeader();
563a673f 375
376 // Get the generator vertex
377 TArrayF mcVertex;
0bd4b00f 378 genHeader->PrimaryVertex(mcVertex);
563a673f 379 Double_t mcVtxZ = mcVertex.At(2);
380
381 // Check the MC vertex is in range
0bd4b00f 382 Int_t mcVtxBin = fVtxAxis.FindBin(mcVtxZ);
383 if (mcVtxBin < 1 || mcVtxBin > fVtxAxis.GetNbins()) {
563a673f 384#ifdef VERBOSE
385 AliWarning(Form("MC v_z=%f is out of rante [%,%f]",
0bd4b00f 386 mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
563a673f 387#endif
388 return;
389 }
390
0bd4b00f 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
563a673f 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);
0bd4b00f 406#endif
407
563a673f 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();
0bd4b00f 434 Double_t phi = particle->Phi();
563a673f 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
0bd4b00f 440 if (!gotInel || gotVertex) continue;
563a673f 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
0bd4b00f 447 if (!trackRef) continue;
563a673f 448
449 // Check that we hit an FMD element
0bd4b00f 450 if (trackRef->DetectorId() != AliTrackReference::kFMD)
563a673f 451 continue;
452
453 // Get the detector coordinates
454 UShort_t d, s, t;
455 Char_t r;
0bd4b00f 456 AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t);
563a673f 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
0bd4b00f 466 // Double_t trRefEta = esd->GetFMDData()->Eta(d,r,s,t);
467 // Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
563a673f 468
469 // Fill strip information into histograms
0bd4b00f 470 FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
563a673f 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//____________________________________________________________________
485void
0bd4b00f 486AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx,
563a673f 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//____________________________________________________________________
502void
0bd4b00f 503AliForwardMCCorrectionsTask::FillStrip(UShort_t d, Char_t r,
563a673f 504 Double_t vz, Double_t eta, Double_t phi,
505 Bool_t first)
506{
507 // Number of hit strips per eta bin
0bd4b00f 508 TH1D* strips = 0;
563a673f 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)
0bd4b00f 512 TH3D* hits = 0;
563a673f 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//____________________________________________________________________
535TH2D*
0bd4b00f 536AliForwardMCCorrectionsTask::GetVertexProj(Int_t v, TH3D* src) const
563a673f 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//____________________________________________________________________
554void
0bd4b00f 555AliForwardMCCorrectionsTask::Terminate(Option_t*)
563a673f 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 =
0bd4b00f 573 static_cast<TH1I*>(list->FindObject(GetEventName(false,false)));
563a673f 574 TH1I* eventsTr =
0bd4b00f 575 static_cast<TH1I*>(list->FindObject(GetEventName(true,false)));
563a673f 576 TH1I* eventsVtx =
0bd4b00f 577 static_cast<TH1I*>(list->FindObject(GetEventName(false,true)));
563a673f 578 TH1I* eventsTrVtx =
0bd4b00f 579 static_cast<TH1I*>(list->FindObject(GetEventName(true,true)));
563a673f 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",
0bd4b00f 628 fVtxAxis.GetNbins(),
629 fVtxAxis.GetXmin(),
630 fVtxAxis.GetXmax());
563a673f 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;
0bd4b00f 643 vl->SetName(Form("vtx%02d", v));
563a673f 644 list->Add(vl);
645
646 // Get event counts
0bd4b00f 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);
563a673f 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
0bd4b00f 662 primIAllV->Scale(1. / nEventsAll);
663 primOAllV->Scale(1. / nEventsAll);
664 primITrVtxV->Scale(1. / nEventsTr);
665 primOTrVtxV->Scale(1. / nEventsTr);
563a673f 666
667 // Calculate the vertex bias on the d^2N/detadphi
0bd4b00f 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)));
563a673f 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)
0bd4b00f 702 TH1D* hStrips = 0;
563a673f 703 switch (d) {
0bd4b00f 704 case 1: hStrips = strips1i; break;
705 case 2: hStrips = (q == 0 ? strips2i : strips2o); break;
706 case 3: hStrips = (q == 0 ? strips3i : strips3o); break;
563a673f 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();
0bd4b00f 716 doubleHit->Divide(hStrips);
563a673f 717 list->Add(doubleHit);
718 }
719 }
720 }
721}
722
723//____________________________________________________________________
724void
0bd4b00f 725AliForwardMCCorrectionsTask::Print(Option_t*) const
563a673f 726{
727}
728
729//
730// EOF
731//