Minor things
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliCentralMultiplicityTask.cxx
CommitLineData
6f791cc3 1//====================================================================
2//
3// Base class for classes that calculate the multiplicity in the
4// central region event-by-event
5//
6// Inputs:
7// - AliESDEvent
8//
9// Outputs:
10// - AliAODCentralMult
11//
12// Histograms
13//
14// Corrections used
15#include "AliCentralMultiplicityTask.h"
52047b6f 16#include "AliAODForwardMult.h"
6f791cc3 17#include "AliForwardUtil.h"
18#include "AliLog.h"
19#include "AliAODHandler.h"
6f791cc3 20#include "AliAnalysisManager.h"
21#include "AliESDEvent.h"
22#include "AliMultiplicity.h"
23#include <TROOT.h>
24#include <TFile.h>
3e478dba 25#include <TError.h>
6f791cc3 26#include <iostream>
27#include <iomanip>
28
29//====================================================================
30AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
31 : AliAnalysisTaskSE(name),
52047b6f 32 fInspector("centralEventInspector"),
6f791cc3 33 fData(0),
34 fList(0),
35 fAODCentral(kFALSE),
3b2bfb07 36 fManager(),
e58000b7 37 fUseSecondary(true),
9453b19e 38 fUseAcceptance(true),
52047b6f 39 fFirstEventSeen(false),
40 fIvz(0)
6f791cc3 41{
fb3430ac 42 //
43 // Constructor
44 //
6f791cc3 45 DefineOutput(1, TList::Class());
a59cbd24 46 fBranchNames =
47 "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
48 "SPDVertex.,PrimaryVertex.";
6f791cc3 49}
50//____________________________________________________________________
9c825779 51AliCentralMultiplicityTask::AliCentralMultiplicityTask()
52 : AliAnalysisTaskSE(),
52047b6f 53 fInspector(),
9c825779 54 fData(0),
55 fList(0),
56 fAODCentral(),
3b2bfb07 57 fManager(),
e58000b7 58 fUseSecondary(true),
9453b19e 59 fUseAcceptance(true),
52047b6f 60 fFirstEventSeen(false),
61 fIvz(0)
9c825779 62{
fb3430ac 63 //
64 // Constructor
65 //
9c825779 66}
67//____________________________________________________________________
68AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
69 : AliAnalysisTaskSE(o),
52047b6f 70 fInspector(o.fInspector),
9c825779 71 fData(o.fData),
72 fList(o.fList),
73 fAODCentral(o.fAODCentral),
3b2bfb07 74 fManager(o.fManager),
e58000b7 75 fUseSecondary(o.fUseSecondary),
9453b19e 76 fUseAcceptance(o.fUseAcceptance),
52047b6f 77 fFirstEventSeen(o.fFirstEventSeen),
78 fIvz(0)
9c825779 79{
fb3430ac 80 //
81 // Copy constructor
82 //
9c825779 83}
84//____________________________________________________________________
85AliCentralMultiplicityTask&
86AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
87{
fb3430ac 88 //
89 // Assignment operator
90 //
52047b6f 91 fInspector = o.fInspector;
9453b19e 92 fData = o.fData;
93 fList = o.fList;
94 fAODCentral = o.fAODCentral;
95 fManager = o.fManager;
96 fUseSecondary = o.fUseSecondary;
97 fUseAcceptance = o.fUseAcceptance;
fb3430ac 98 fFirstEventSeen = o.fFirstEventSeen;
52047b6f 99 fIvz = 0;
9c825779 100 return *this;
101}
102//____________________________________________________________________
3e478dba 103void AliCentralMultiplicityTask::UserCreateOutputObjects()
104{
fb3430ac 105 //
106 // Create output objects
107 //
108 //
6f791cc3 109
110 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
111 AliAODHandler* ah =
112 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
113 if (!ah) AliFatal("No AOD output handler set in analysis manager");
114
115
116 TObject* obj = &fAODCentral;
117 ah->AddBranch("AliAODCentralMult", &obj);
52047b6f 118
119
6f791cc3 120 fList = new TList();
9d05ffeb 121 fList->SetOwner();
52047b6f 122
123 fInspector.DefineOutput(fList);
124
125 PostData(1,fList);
126}
127
128//____________________________________________________________________
129AliESDEvent*
130AliCentralMultiplicityTask::GetESDEvent()
131{
132 //
133 // Get the ESD event. IF this is the first event, initialise
134 //
135 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
136 if (!esd) {
137 AliWarning("No ESD event found for input event");
138 return 0;
139 }
140
141 if (fFirstEventSeen) return esd;
142 if (GetManager().IsInit()) return esd;
6f791cc3 143
52047b6f 144 fInspector.ReadRunDetails(esd);
145 GetManager().Init(fInspector.GetCollisionSystem(),
146 fInspector.GetEnergy(),
147 fInspector.GetField());
148 AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
149 if (!secMap)
150 AliFatal("No secondary map defined!");
151 const TAxis& vaxis = secMap->GetVertexAxis();
152
153 std::cout << "Vertex range is "
154 << vaxis.GetNbins() << ","
155 << vaxis.GetXmin() << ","
156 << vaxis.GetXmax()
157 << std::endl;
158 fInspector.Init(vaxis);
159 AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
160 fFirstEventSeen = kTRUE;
161 Print();
162
163 return esd;
6f791cc3 164}
52047b6f 165//____________________________________________________________________
166void
167AliCentralMultiplicityTask::MarkEventForStore() const
168{
169 // Make sure the AOD tree is filled
170 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
171 AliAODHandler* ah =
172 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
173 if (!ah)
174 AliFatal("No AOD output handler set in analysis manager");
175
176 ah->SetFillAOD(kTRUE);
177}
178
6f791cc3 179//____________________________________________________________________
3e478dba 180void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
181{
fb3430ac 182 //
183 // Process each event
184 //
185 // Parameters:
186 // option Not used
187 //
c4a7e081 188 fAODCentral.Clear("");
52047b6f 189 fIvz = 0;
190
191 AliESDEvent* esd = GetESDEvent();
6f791cc3 192
52047b6f 193 Bool_t lowFlux = kFALSE;
194 UInt_t triggers = 0;
195 UShort_t ivz = 0;
196 Double_t vz = 0;
197 Double_t cent = -1;
198 UShort_t nClusters = 0;
199 UInt_t found = fInspector.Process(esd, triggers, lowFlux,
200 ivz, vz, cent, nClusters);
201
202 // No event or no trigger
203 if (found & AliFMDEventInspector::kNoEvent) return;
204 if (found & AliFMDEventInspector::kNoTriggers) return;
6f791cc3 205
206 // Make sure AOD is filled
52047b6f 207 MarkEventForStore();
208
209 if (found == AliFMDEventInspector::kNoSPD) return;
210 if (found == AliFMDEventInspector::kNoVertex) return;
211 if (triggers & AliAODForwardMult::kPileUp) return;
212 if (found == AliFMDEventInspector::kBadVertex) return; // Out of range
6f791cc3 213
214 //Doing analysis
52047b6f 215 fIvz = ivz;
6f791cc3 216 const AliMultiplicity* spdmult = esd->GetMultiplicity();
52047b6f 217
218 TH2D& aodHist = fAODCentral.GetHistogram();
219
220 ProcessESD(aodHist, spdmult);
221 CorrectData(aodHist, ivz);
222
223 PostData(1,fList);
224}
225//____________________________________________________________________
226void
227AliCentralMultiplicityTask::ProcessESD(TH2D& aodHist,
228 const AliMultiplicity* spdmult) const
229{
230
6f791cc3 231 //Filling clusters in layer 1 used for tracklets...
3e478dba 232 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
52047b6f 233 aodHist.Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
3e478dba 234
6f791cc3 235 //...and then the unused ones in layer 1
3e478dba 236 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
52047b6f 237 aodHist.Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
238 spdmult->GetPhiSingle(j));
239}
240
241//____________________________________________________________________
242void
243AliCentralMultiplicityTask::CorrectData(TH2D& aodHist, UShort_t vtxbin) const
244{
3e478dba 245 // Corrections
6f791cc3 246 TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
9453b19e 247 TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
12fffad7 248
9453b19e 249 if (!hSecMap) AliFatal("No secondary map!");
250 if (!hAcceptance) AliFatal("No acceptance!");
dbe8d9ed 251
52047b6f 252 if (fUseSecondary && hSecMap) aodHist.Divide(hSecMap);
6f791cc3 253
52047b6f 254 for(Int_t nx = 1; nx <= aodHist.GetNbinsX(); nx++) {
3b2bfb07 255 Float_t accCor = hAcceptance->GetBinContent(nx);
9453b19e 256 Float_t accErr = hAcceptance->GetBinError(nx);
6f791cc3 257
258 Bool_t etabinSeen = kFALSE;
52047b6f 259 for(Int_t ny = 1; ny <= aodHist.GetNbinsY(); ny++) {
9453b19e 260 // Get currrent value
52047b6f 261 Float_t aodValue = aodHist.GetBinContent(nx,ny);
262 Float_t aodErr = aodHist.GetBinError(nx,ny);
9453b19e 263
264 // Set underflow bin
12fffad7 265 Float_t secCor = 0;
52047b6f 266 if(hSecMap) secCor = hSecMap->GetBinContent(nx,ny);
3b2bfb07 267 if (secCor > 0.5) etabinSeen = kTRUE;
52047b6f 268 if (aodValue < 0.000001) {
269 aodHist.SetBinContent(nx,ny, 0);
270 continue;
271 }
9453b19e 272
273 if (!fUseAcceptance) continue;
274
275 // Acceptance correction
3b2bfb07 276 if (accCor < 0.000001) accCor = 1;
277 Float_t aodNew = aodValue / accCor ;
12fffad7 278 Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
279 TMath::Power(accErr/accCor,2) );
52047b6f 280 aodHist.SetBinContent(nx,ny, aodNew);
12fffad7 281 //test
52047b6f 282 aodHist.SetBinError(nx,ny,error);
283 aodHist.SetBinError(nx,ny,aodErr);
6f791cc3 284
285 }
286 //Filling underflow bin if we eta bin is in range
52047b6f 287 if(etabinSeen) aodHist.SetBinContent(nx,0, 1.);
6f791cc3 288 }
6f791cc3 289}
52047b6f 290
6f791cc3 291//____________________________________________________________________
3e478dba 292void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
293{
fb3430ac 294 //
295 // End of job
296 //
297 // Parameters:
298 // option Not used
299 //
6f791cc3 300}
301//____________________________________________________________________
302void
52047b6f 303AliCentralMultiplicityTask::Print(Option_t* option) const
6f791cc3 304{
fb3430ac 305 //
306 // Print information
307 //
308 // Parameters:
309 // option Not used
310 //
52047b6f 311 std::cout << ClassName() << ": " << GetName() << "\n"
312 << std::boolalpha
313 << " Use secondary correction: " << fUseSecondary << '\n'
314 << " Use acceptance correction: " << fUseAcceptance << '\n'
315 << " Off-line trigger mask: 0x"
316 << std::hex << std::setfill('0')
317 << std::setw (8) << fOfflineTriggerMask
318 << std::dec << std::setfill (' ')
319 << std::noboolalpha << std::endl;
320 gROOT->IncreaseDirLevel();
321 fManager.Print(option);
322 fInspector.Print(option);
323 gROOT->DecreaseDirLevel();
324
6f791cc3 325}
3e478dba 326//====================================================================
6f791cc3 327AliCentralMultiplicityTask::Manager::Manager() :
328 fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
329 fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
330 fAcceptance(),
331 fSecmap(),
332 fAcceptanceName("centralacceptance"),
e58000b7 333 fSecMapName("centralsecmap"),
334 fIsInit(kFALSE)
6f791cc3 335{
fb3430ac 336 //
337 // Constructor
338 //
6f791cc3 339}
340//____________________________________________________________________
3e478dba 341AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
342 :fAcceptancePath(o.fAcceptancePath),
343 fSecMapPath(o.fSecMapPath),
344 fAcceptance(o.fAcceptance),
345 fSecmap(o.fSecmap),
346 fAcceptanceName(o.fAcceptanceName),
e58000b7 347 fSecMapName(o.fSecMapName),
348 fIsInit(o.fIsInit)
fb3430ac 349{
350 //
351 // Copy Constructor
352 //
353}
3e478dba 354//____________________________________________________________________
355AliCentralMultiplicityTask::Manager&
356AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
357{
fb3430ac 358 //
359 // Assignment operator
360 //
3e478dba 361 fAcceptancePath = o.fAcceptancePath;
362 fSecMapPath = o.fSecMapPath;
363 fAcceptance = o.fAcceptance;
364 fSecmap = o.fSecmap;
365 fAcceptanceName = o.fAcceptanceName;
366 fSecMapName = o.fSecMapName;
e58000b7 367 fIsInit = o.fIsInit;
3e478dba 368 return *this;
369}
370
371//____________________________________________________________________
372const char*
373AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
374 UShort_t sys,
375 UShort_t sNN,
376 Short_t field) const
377{
fb3430ac 378 //
379 // Get full path name to object file
380 //
381 // Parameters:
382 // what What to get
383 // sys Collision system
384 // sNN Center of mass energy
385 // field Magnetic field
386 //
387 // Return:
388 //
389 //
3e478dba 390 return Form("%s/%s",
391 what == 0 ? GetSecMapPath() : GetAcceptancePath(),
392 GetFileName(what, sys, sNN, field));
393}
394
395//____________________________________________________________________
396const char*
397AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
398 UShort_t sys,
399 UShort_t sNN,
400 Short_t field) const
401{
fb3430ac 402 //
403 // Get the full path name
404 //
405 // Parameters:
406 // what What to get
407 // sys Collision system
408 // sNN Center of mass energy
409 // field Magnetic field
410 //
411 // Return:
412 //
413 //
3e478dba 414 // Must be static - otherwise the data may disappear on return from
415 // this member function
416 static TString fname = "";
6f791cc3 417
dbe8d9ed 418 fname = "";
6f791cc3 419 switch(what) {
3e478dba 420 case 0: fname.Append(fSecMapName.Data()); break;
421 case 1: fname.Append(fAcceptanceName.Data()); break;
6f791cc3 422 default:
3e478dba 423 ::Error("GetFileName",
424 "Invalid indentifier %d for central object, must be 0 or 1!", what);
6f791cc3 425 break;
426 }
427 fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
428 AliForwardUtil::CollisionSystemString(sys),
429 sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
430
431 return fname.Data();
6f791cc3 432}
3e478dba 433
6f791cc3 434//____________________________________________________________________
3e478dba 435TH2D*
436AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
437{
fb3430ac 438 //
439 // Get the secondary map
440 //
441 // Parameters:
442 // vtxbin
443 //
444 // Return:
445 //
446 //
3e478dba 447 if (!fSecmap) {
448 ::Warning("GetSecMapCorrection","No secondary map defined");
449 return 0;
450 }
451 return fSecmap->GetCorrection(vtxbin);
452}
453//____________________________________________________________________
454TH1D*
455AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
456 const
457{
fb3430ac 458 //
459 // Get the acceptance correction
460 //
461 // Parameters:
462 // vtxbin
463 //
464 // Return:
465 //
466 //
3e478dba 467 if (!fAcceptance) {
468 ::Warning("GetAcceptanceCorrection","No acceptance map defined");
469 return 0;
470 }
471 return fAcceptance->GetCorrection(vtxbin);
472}
473
474//____________________________________________________________________
475void
476AliCentralMultiplicityTask::Manager::Init(UShort_t sys,
477 UShort_t sNN,
478 Short_t field)
479{
fb3430ac 480 //
481 // Initialize
482 //
483 // Parameters:
484 // sys Collision system (1: pp, 2: PbPb)
485 // sNN Center of mass energy per nucleon pair [GeV]
486 // field Magnetic field [kG]
487 //
e58000b7 488 if(fIsInit) ::Warning("Init","Already initialised - overriding...");
489
6f791cc3 490 TFile fsec(GetFullFileName(0,sys,sNN,field));
3e478dba 491 fSecmap =
492 dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
6f791cc3 493 if(!fSecmap) {
3e478dba 494 ::Error("Init", "no central Secondary Map found!") ;
6f791cc3 495 return;
496 }
497 TFile facc(GetFullFileName(1,sys,sNN,field));
3e478dba 498 fAcceptance =
499 dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
500 if(!fAcceptance) {
501 ::Error("Init", "no central Acceptance found!") ;
502 return;
503 }
e58000b7 504
505 if(fSecmap && fAcceptance) {
506 fIsInit = kTRUE;
507 ::Info("Init",
52047b6f 508 "Central Manager initialised for %s, energy %dGeV, field %dkG",
509 sys == 1 ? "pp" : sys == 2 ? "PbPb" : "unknown", sNN,field);
510 }
511}
6f791cc3 512
52047b6f 513//____________________________________________________________________
514void
515AliCentralMultiplicityTask::Manager::Print(Option_t* option) const
516{
517 std::cout << " AliCentralMultiplicityTask::Manager\n"
518 << std::boolalpha
519 << " Initialized: " << fIsInit << '\n'
520 << " Acceptance path: " << fAcceptancePath << '\n'
521 << " Acceptance name: " << fAcceptanceName << '\n'
522 << " Acceptance: " << fAcceptance << '\n'
523 << " Secondary path: " << fSecMapPath << '\n'
524 << " Secondary name: " << fSecMapName << '\n'
525 << " Secondary map: " << fSecmap
526 << std::noboolalpha << std::endl;
527 if (fAcceptance) fAcceptance->Print(option);
528 if (fSecmap) fSecmap->Print(option);
6f791cc3 529}
52047b6f 530
6f791cc3 531//
532// EOF
533//