]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Various improvements
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMCCorrectionsTask.cxx
CommitLineData
7984e5f7 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//
0bd4b00f 14#include "AliForwardMCCorrectionsTask.h"
563a673f 15#include "AliTriggerAnalysis.h"
16#include "AliPhysicsSelection.h"
17#include "AliLog.h"
0bd4b00f 18#include "AliHeader.h"
19#include "AliGenEventHeader.h"
563a673f 20#include "AliESDEvent.h"
21#include "AliAODHandler.h"
22#include "AliMultiplicity.h"
23#include "AliInputEventHandler.h"
0bd4b00f 24#include "AliStack.h"
25#include "AliMCEvent.h"
26#include "AliFMDStripIndex.h"
563a673f 27#include <TH1.h>
0bd4b00f 28#include <TH3D.h>
563a673f 29#include <TDirectory.h>
30#include <TTree.h>
31
32//====================================================================
33namespace {
34 const char* GetEventName(Bool_t tr, Bool_t vtx)
35 {
0bd4b00f 36 return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
563a673f 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//====================================================================
0bd4b00f 55AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
563a673f 56 : AliAnalysisTaskSE(),
0bd4b00f 57 fHEvents(0),
563a673f 58 fHEventsTr(0),
59 fHEventsTrVtx(0),
0bd4b00f 60 fHEventsVtx(0),
563a673f 61 fHTriggers(0),
0bd4b00f 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),
563a673f 76 fVtxAxis(),
0bd4b00f 77 fEtaAxis(),
78 fList()
563a673f 79{
7984e5f7 80 //
81 // Constructor
82 //
83 // Parameters:
84 // name Name of task
85 //
563a673f 86}
87
88//____________________________________________________________________
0bd4b00f 89AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
563a673f 90 : AliAnalysisTaskSE(name),
0bd4b00f 91 fHEvents(0),
563a673f 92 fHEventsTr(0),
0bd4b00f 93 fHEventsTrVtx(0),
94 fHEventsVtx(0),
563a673f 95 fHTriggers(0),
0bd4b00f 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),
563a673f 110 fVtxAxis(10,-10,10),
0bd4b00f 111 fEtaAxis(200,-4,6),
112 fList()
563a673f 113{
7984e5f7 114 //
115 // Constructor
116 //
117 // Parameters:
118 // name Name of task
119 //
563a673f 120 DefineOutput(1, TList::Class());
121}
122
123//____________________________________________________________________
0bd4b00f 124AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
563a673f 125 : AliAnalysisTaskSE(o),
0bd4b00f 126 fHEvents(o.fHEvents),
563a673f 127 fHEventsTr(o.fHEventsTr),
0bd4b00f 128 fHEventsTrVtx(o.fHEventsTrVtx),
129 fHEventsVtx(o.fHEventsVtx),
563a673f 130 fHTriggers(o.fHTriggers),
0bd4b00f 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)
563a673f 148{
7984e5f7 149 //
150 // Copy constructor
151 //
152 // Parameters:
153 // o Object to copy from
154 //
563a673f 155}
156
157//____________________________________________________________________
0bd4b00f 158AliForwardMCCorrectionsTask&
159AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
563a673f 160{
7984e5f7 161 //
162 // Assignment operator
163 //
164 // Parameters:
165 // o Object to assign from
166 //
167 // Return:
168 // Reference to this object
169 //
563a673f 170 fHEventsTr = o.fHEventsTr;
171 fHEventsTrVtx = o.fHEventsTrVtx;
172 fHTriggers = o.fHTriggers;
0bd4b00f 173 SetVertexAxis(o.fVtxAxis);
174 SetEtaAxis(o.fEtaAxis);
563a673f 175
176 return *this;
177}
178
179//____________________________________________________________________
180void
0bd4b00f 181AliForwardMCCorrectionsTask::Init()
563a673f 182{
7984e5f7 183 //
184 // Initialize the task
185 //
186 //
563a673f 187}
188
189//____________________________________________________________________
190void
0bd4b00f 191AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
192 Double_t max)
563a673f 193{
7984e5f7 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 //
563a673f 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//____________________________________________________________________
215void
0bd4b00f 216AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
563a673f 217{
7984e5f7 218 //
219 // Set the vertex axis to use
220 //
221 // Parameters:
222 // axis Axis
223 //
563a673f 224 SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
225}
226
227//____________________________________________________________________
228void
0bd4b00f 229AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
563a673f 230{
7984e5f7 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 //
563a673f 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//____________________________________________________________________
252void
0bd4b00f 253AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
563a673f 254{
7984e5f7 255 //
256 // Set the eta axis to use
257 //
258 // Parameters:
259 // axis Axis
260 //
563a673f 261 SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
262}
263
563a673f 264//____________________________________________________________________
265TH3D*
0bd4b00f 266AliForwardMCCorrectionsTask::Make3D(const char* name, const char* title,
563a673f 267 Int_t nPhi) const
268{
7984e5f7 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 //
563a673f 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//____________________________________________________________________
298TH1D*
0bd4b00f 299AliForwardMCCorrectionsTask::Make1D(const char* name, const char* title) const
563a673f 300{
7984e5f7 301 //
302 // Make 1D histogram
303 //
304 // Parameters:
305 // name Name
306 // title Title
307 //
308 // Return:
309 // Histogram
310 //
563a673f 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//____________________________________________________________________
326void
0bd4b00f 327AliForwardMCCorrectionsTask::UserCreateOutputObjects()
563a673f 328{
7984e5f7 329 //
330 // Create output objects
331 //
332 //
563a673f 333 fList = new TList;
334 fList->SetName(GetName());
335
0bd4b00f 336 fHEvents = new TH1I(GetEventName(false,false),
563a673f 337 "Number of all events",
0bd4b00f 338 fVtxAxis.GetNbins(),
339 fVtxAxis.GetXmin(),
340 fVtxAxis.GetXmax());
563a673f 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
0bd4b00f 348 fHEventsTr = new TH1I(GetEventName(true, false),
349 "Number of triggered events",
563a673f 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
0bd4b00f 360 fHEventsTrVtx = new TH1I(GetEventName(true,true),
563a673f 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);
0bd4b00f 371
372 fHEventsVtx = new TH1I(GetEventName(false,true),
563a673f 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//____________________________________________________________________
448void
0bd4b00f 449AliForwardMCCorrectionsTask::UserExec(Option_t*)
563a673f 450{
7984e5f7 451 //
452 // Process each event
453 //
454 // Parameters:
455 // option Not used
456 //
457
563a673f 458 // Get the input data - MC event
459 AliMCEvent* mcEvent = MCEvent();
cc83fca2 460 if (!mcEvent) {
563a673f 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();
0bd4b00f 477 AliGenEventHeader* genHeader = mcHeader->GenEventHeader();
563a673f 478
479 // Get the generator vertex
480 TArrayF mcVertex;
0bd4b00f 481 genHeader->PrimaryVertex(mcVertex);
563a673f 482 Double_t mcVtxZ = mcVertex.At(2);
483
484 // Check the MC vertex is in range
0bd4b00f 485 Int_t mcVtxBin = fVtxAxis.FindBin(mcVtxZ);
486 if (mcVtxBin < 1 || mcVtxBin > fVtxAxis.GetNbins()) {
563a673f 487#ifdef VERBOSE
488 AliWarning(Form("MC v_z=%f is out of rante [%,%f]",
0bd4b00f 489 mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
563a673f 490#endif
491 return;
492 }
493
0bd4b00f 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
563a673f 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);
0bd4b00f 509#endif
510
563a673f 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();
0bd4b00f 537 Double_t phi = particle->Phi();
563a673f 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
0bd4b00f 543 if (!gotInel || gotVertex) continue;
563a673f 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
0bd4b00f 550 if (!trackRef) continue;
563a673f 551
552 // Check that we hit an FMD element
0bd4b00f 553 if (trackRef->DetectorId() != AliTrackReference::kFMD)
563a673f 554 continue;
555
556 // Get the detector coordinates
557 UShort_t d, s, t;
558 Char_t r;
0bd4b00f 559 AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t);
563a673f 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
0bd4b00f 569 // Double_t trRefEta = esd->GetFMDData()->Eta(d,r,s,t);
570 // Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
563a673f 571
572 // Fill strip information into histograms
0bd4b00f 573 FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
563a673f 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//____________________________________________________________________
588void
0bd4b00f 589AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx,
563a673f 590 Double_t vz, Double_t eta, Double_t phi)
591{
7984e5f7 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 //
563a673f 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//____________________________________________________________________
615void
0bd4b00f 616AliForwardMCCorrectionsTask::FillStrip(UShort_t d, Char_t r,
563a673f 617 Double_t vz, Double_t eta, Double_t phi,
618 Bool_t first)
619{
7984e5f7 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
563a673f 632 // Number of hit strips per eta bin
0bd4b00f 633 TH1D* strips = 0;
563a673f 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)
0bd4b00f 637 TH3D* hits = 0;
563a673f 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//____________________________________________________________________
660TH2D*
0bd4b00f 661AliForwardMCCorrectionsTask::GetVertexProj(Int_t v, TH3D* src) const
563a673f 662{
7984e5f7 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 //
563a673f 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//____________________________________________________________________
689void
0bd4b00f 690AliForwardMCCorrectionsTask::Terminate(Option_t*)
563a673f 691{
7984e5f7 692 //
693 // End of job
694 //
695 // Parameters:
696 // option Not used
697 //
698
563a673f 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 =
0bd4b00f 715 static_cast<TH1I*>(list->FindObject(GetEventName(false,false)));
563a673f 716 TH1I* eventsTr =
0bd4b00f 717 static_cast<TH1I*>(list->FindObject(GetEventName(true,false)));
563a673f 718 TH1I* eventsVtx =
0bd4b00f 719 static_cast<TH1I*>(list->FindObject(GetEventName(false,true)));
563a673f 720 TH1I* eventsTrVtx =
0bd4b00f 721 static_cast<TH1I*>(list->FindObject(GetEventName(true,true)));
563a673f 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",
0bd4b00f 770 fVtxAxis.GetNbins(),
771 fVtxAxis.GetXmin(),
772 fVtxAxis.GetXmax());
563a673f 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;
0bd4b00f 785 vl->SetName(Form("vtx%02d", v));
563a673f 786 list->Add(vl);
787
788 // Get event counts
0bd4b00f 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);
563a673f 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
0bd4b00f 804 primIAllV->Scale(1. / nEventsAll);
805 primOAllV->Scale(1. / nEventsAll);
806 primITrVtxV->Scale(1. / nEventsTr);
807 primOTrVtxV->Scale(1. / nEventsTr);
563a673f 808
809 // Calculate the vertex bias on the d^2N/detadphi
0bd4b00f 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)));
563a673f 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)
0bd4b00f 844 TH1D* hStrips = 0;
563a673f 845 switch (d) {
0bd4b00f 846 case 1: hStrips = strips1i; break;
847 case 2: hStrips = (q == 0 ? strips2i : strips2o); break;
848 case 3: hStrips = (q == 0 ? strips3i : strips3o); break;
563a673f 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();
0bd4b00f 858 doubleHit->Divide(hStrips);
563a673f 859 list->Add(doubleHit);
860 }
861 }
862 }
863}
864
865//____________________________________________________________________
866void
0bd4b00f 867AliForwardMCCorrectionsTask::Print(Option_t*) const
563a673f 868{
869}
870
871//
872// EOF
873//