Bug fix in vertex recalculation
[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"
16#include "AliForwardCorrectionManager.h"
17#include "AliForwardUtil.h"
18#include "AliLog.h"
19#include "AliAODHandler.h"
20#include "AliInputEventHandler.h"
21#include "AliESDInputHandler.h"
22#include "AliAnalysisManager.h"
23#include "AliESDEvent.h"
24#include "AliMultiplicity.h"
e58000b7 25#include "AliFMDEventInspector.h"
6f791cc3 26#include <TROOT.h>
27#include <TFile.h>
3e478dba 28#include <TError.h>
6f791cc3 29#include <iostream>
30#include <iomanip>
31
32//====================================================================
33AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
34 : AliAnalysisTaskSE(name),
35 fData(0),
36 fList(0),
37 fAODCentral(kFALSE),
3b2bfb07 38 fManager(),
e58000b7 39 fUseSecondary(true),
9453b19e 40 fUseAcceptance(true),
fb3430ac 41 fFirstEventSeen(false)
6f791cc3 42{
fb3430ac 43 //
44 // Constructor
45 //
6f791cc3 46 DefineOutput(1, TList::Class());
a59cbd24 47 fBranchNames =
48 "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
49 "SPDVertex.,PrimaryVertex.";
6f791cc3 50}
51//____________________________________________________________________
9c825779 52AliCentralMultiplicityTask::AliCentralMultiplicityTask()
53 : AliAnalysisTaskSE(),
54 fData(0),
55 fList(0),
56 fAODCentral(),
3b2bfb07 57 fManager(),
e58000b7 58 fUseSecondary(true),
9453b19e 59 fUseAcceptance(true),
fb3430ac 60 fFirstEventSeen(false)
9c825779 61{
fb3430ac 62 //
63 // Constructor
64 //
9c825779 65}
66//____________________________________________________________________
67AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
68 : AliAnalysisTaskSE(o),
69 fData(o.fData),
70 fList(o.fList),
71 fAODCentral(o.fAODCentral),
3b2bfb07 72 fManager(o.fManager),
e58000b7 73 fUseSecondary(o.fUseSecondary),
9453b19e 74 fUseAcceptance(o.fUseAcceptance),
fb3430ac 75 fFirstEventSeen(o.fFirstEventSeen)
9c825779 76{
fb3430ac 77 //
78 // Copy constructor
79 //
9c825779 80}
81//____________________________________________________________________
82AliCentralMultiplicityTask&
83AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
84{
fb3430ac 85 //
86 // Assignment operator
87 //
9453b19e 88 fData = o.fData;
89 fList = o.fList;
90 fAODCentral = o.fAODCentral;
91 fManager = o.fManager;
92 fUseSecondary = o.fUseSecondary;
93 fUseAcceptance = o.fUseAcceptance;
fb3430ac 94 fFirstEventSeen = o.fFirstEventSeen;
9c825779 95 return *this;
96}
97//____________________________________________________________________
3e478dba 98void AliCentralMultiplicityTask::UserCreateOutputObjects()
99{
fb3430ac 100 //
101 // Create output objects
102 //
103 //
6f791cc3 104
105 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
106 AliAODHandler* ah =
107 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
108 if (!ah) AliFatal("No AOD output handler set in analysis manager");
109
110
111 TObject* obj = &fAODCentral;
112 ah->AddBranch("AliAODCentralMult", &obj);
113
114 fList = new TList();
9d05ffeb 115 fList->SetOwner();
6f791cc3 116 PostData(1,fList);
117
118}
119//____________________________________________________________________
3e478dba 120void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
121{
fb3430ac 122 //
123 // Process each event
124 //
125 // Parameters:
126 // option Not used
127 //
c4a7e081 128 fAODCentral.Clear("");
6f791cc3 129
130 AliESDInputHandler* eventHandler =
131 dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
132 ->GetInputEventHandler());
133 if (!eventHandler) {
134 AliWarning("No inputhandler found for this event!");
135 return;}
136
137 AliESDEvent* esd = eventHandler->GetEvent();
138
fb3430ac 139 if(!GetManager().IsInit() && !fFirstEventSeen) {
e58000b7 140 AliFMDEventInspector inspector;
141 inspector.ReadRunDetails(esd);
142 GetManager().Init(inspector.GetCollisionSystem(),
143 inspector.GetEnergy(),
144 inspector.GetField());
145
e58000b7 146 AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
fb3430ac 147 fFirstEventSeen = kTRUE;
e58000b7 148 }
c4a7e081 149
6f791cc3 150 //Selecting only events with |valid vertex| < 10 cm
151 const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
152 if (!vertex) return;
153 if(!(vertex->GetStatus())) return;
154 if(vertex->GetNContributors() <= 0) return ;
155 if(vertex->GetZRes() > 0.1 ) return;
156 Double_t vertexXYZ[3]={0,0,0};
157 vertex->GetXYZ(vertexXYZ);
158 if(TMath::Abs(vertexXYZ[2]) > 10) return;
159
3e478dba 160 Double_t delta = 2 ;
6f791cc3 161 Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
3e478dba 162 //HHD: The vtxbins are 1-10, not 0-9
163 Int_t vtxbin = Int_t(vertexBinDouble + 1) ;
6f791cc3 164
165 // Make sure AOD is filled
166 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
167 AliAODHandler* ah =
168 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
169 if (!ah)
170 AliFatal("No AOD output handler set in analysis manager");
171
172 ah->SetFillAOD(kTRUE);
173
174 //Doing analysis
c4a7e081 175
6f791cc3 176 TH2D *aodHist = &(fAODCentral.GetHistogram());
177
178 const AliMultiplicity* spdmult = esd->GetMultiplicity();
179 //Filling clusters in layer 1 used for tracklets...
3e478dba 180 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
6f791cc3 181 aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
3e478dba 182
6f791cc3 183 //...and then the unused ones in layer 1
3e478dba 184 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
6f791cc3 185 aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
3e478dba 186 spdmult->GetPhiSingle(j));
6f791cc3 187
6f791cc3 188
3e478dba 189 // Corrections
6f791cc3 190 TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
9453b19e 191 TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
12fffad7 192
9453b19e 193 if (!hSecMap) AliFatal("No secondary map!");
194 if (!hAcceptance) AliFatal("No acceptance!");
dbe8d9ed 195
9453b19e 196 if (fUseSecondary && hSecMap) aodHist->Divide(hSecMap);
6f791cc3 197
198 for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
3b2bfb07 199 Float_t accCor = hAcceptance->GetBinContent(nx);
9453b19e 200 Float_t accErr = hAcceptance->GetBinError(nx);
6f791cc3 201
202 Bool_t etabinSeen = kFALSE;
203 for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
9453b19e 204 // Get currrent value
3b2bfb07 205 Float_t aodValue = aodHist->GetBinContent(nx,ny);
9453b19e 206 Float_t aodErr = aodHist->GetBinError(nx,ny);
207
208 // Set underflow bin
12fffad7 209 Float_t secCor = 0;
210 if(hSecMap) secCor = hSecMap->GetBinContent(nx,ny);
3b2bfb07 211 if (secCor > 0.5) etabinSeen = kTRUE;
212 if (aodValue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
9453b19e 213
214 if (!fUseAcceptance) continue;
215
216 // Acceptance correction
3b2bfb07 217 if (accCor < 0.000001) accCor = 1;
218 Float_t aodNew = aodValue / accCor ;
12fffad7 219 Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
220 TMath::Power(accErr/accCor,2) );
9453b19e 221 aodHist->SetBinContent(nx,ny, aodNew);
12fffad7 222 //test
6f791cc3 223 aodHist->SetBinError(nx,ny,error);
12fffad7 224 aodHist->SetBinError(nx,ny,aodErr);
6f791cc3 225
226 }
227 //Filling underflow bin if we eta bin is in range
228 if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
6f791cc3 229 }
230
231 PostData(1,fList);
6f791cc3 232}
233//____________________________________________________________________
3e478dba 234void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
235{
fb3430ac 236 //
237 // End of job
238 //
239 // Parameters:
240 // option Not used
241 //
6f791cc3 242}
243//____________________________________________________________________
244void
245AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
246{
fb3430ac 247 //
248 // Print information
249 //
250 // Parameters:
251 // option Not used
252 //
6f791cc3 253}
3e478dba 254//====================================================================
6f791cc3 255AliCentralMultiplicityTask::Manager::Manager() :
256 fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
257 fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
258 fAcceptance(),
259 fSecmap(),
260 fAcceptanceName("centralacceptance"),
e58000b7 261 fSecMapName("centralsecmap"),
262 fIsInit(kFALSE)
6f791cc3 263{
fb3430ac 264 //
265 // Constructor
266 //
6f791cc3 267}
268//____________________________________________________________________
3e478dba 269AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
270 :fAcceptancePath(o.fAcceptancePath),
271 fSecMapPath(o.fSecMapPath),
272 fAcceptance(o.fAcceptance),
273 fSecmap(o.fSecmap),
274 fAcceptanceName(o.fAcceptanceName),
e58000b7 275 fSecMapName(o.fSecMapName),
276 fIsInit(o.fIsInit)
fb3430ac 277{
278 //
279 // Copy Constructor
280 //
281}
3e478dba 282//____________________________________________________________________
283AliCentralMultiplicityTask::Manager&
284AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
285{
fb3430ac 286 //
287 // Assignment operator
288 //
3e478dba 289 fAcceptancePath = o.fAcceptancePath;
290 fSecMapPath = o.fSecMapPath;
291 fAcceptance = o.fAcceptance;
292 fSecmap = o.fSecmap;
293 fAcceptanceName = o.fAcceptanceName;
294 fSecMapName = o.fSecMapName;
e58000b7 295 fIsInit = o.fIsInit;
3e478dba 296 return *this;
297}
298
299//____________________________________________________________________
300const char*
301AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
302 UShort_t sys,
303 UShort_t sNN,
304 Short_t field) const
305{
fb3430ac 306 //
307 // Get full path name to object file
308 //
309 // Parameters:
310 // what What to get
311 // sys Collision system
312 // sNN Center of mass energy
313 // field Magnetic field
314 //
315 // Return:
316 //
317 //
3e478dba 318 return Form("%s/%s",
319 what == 0 ? GetSecMapPath() : GetAcceptancePath(),
320 GetFileName(what, sys, sNN, field));
321}
322
323//____________________________________________________________________
324const char*
325AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
326 UShort_t sys,
327 UShort_t sNN,
328 Short_t field) const
329{
fb3430ac 330 //
331 // Get the full path name
332 //
333 // Parameters:
334 // what What to get
335 // sys Collision system
336 // sNN Center of mass energy
337 // field Magnetic field
338 //
339 // Return:
340 //
341 //
3e478dba 342 // Must be static - otherwise the data may disappear on return from
343 // this member function
344 static TString fname = "";
6f791cc3 345
dbe8d9ed 346 fname = "";
6f791cc3 347 switch(what) {
3e478dba 348 case 0: fname.Append(fSecMapName.Data()); break;
349 case 1: fname.Append(fAcceptanceName.Data()); break;
6f791cc3 350 default:
3e478dba 351 ::Error("GetFileName",
352 "Invalid indentifier %d for central object, must be 0 or 1!", what);
6f791cc3 353 break;
354 }
355 fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
356 AliForwardUtil::CollisionSystemString(sys),
357 sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
358
359 return fname.Data();
6f791cc3 360}
3e478dba 361
6f791cc3 362//____________________________________________________________________
3e478dba 363TH2D*
364AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
365{
fb3430ac 366 //
367 // Get the secondary map
368 //
369 // Parameters:
370 // vtxbin
371 //
372 // Return:
373 //
374 //
3e478dba 375 if (!fSecmap) {
376 ::Warning("GetSecMapCorrection","No secondary map defined");
377 return 0;
378 }
379 return fSecmap->GetCorrection(vtxbin);
380}
381//____________________________________________________________________
382TH1D*
383AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
384 const
385{
fb3430ac 386 //
387 // Get the acceptance correction
388 //
389 // Parameters:
390 // vtxbin
391 //
392 // Return:
393 //
394 //
3e478dba 395 if (!fAcceptance) {
396 ::Warning("GetAcceptanceCorrection","No acceptance map defined");
397 return 0;
398 }
399 return fAcceptance->GetCorrection(vtxbin);
400}
401
402//____________________________________________________________________
403void
404AliCentralMultiplicityTask::Manager::Init(UShort_t sys,
405 UShort_t sNN,
406 Short_t field)
407{
fb3430ac 408 //
409 // Initialize
410 //
411 // Parameters:
412 // sys Collision system (1: pp, 2: PbPb)
413 // sNN Center of mass energy per nucleon pair [GeV]
414 // field Magnetic field [kG]
415 //
e58000b7 416 if(fIsInit) ::Warning("Init","Already initialised - overriding...");
417
6f791cc3 418 TFile fsec(GetFullFileName(0,sys,sNN,field));
3e478dba 419 fSecmap =
420 dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
6f791cc3 421 if(!fSecmap) {
3e478dba 422 ::Error("Init", "no central Secondary Map found!") ;
6f791cc3 423 return;
424 }
425 TFile facc(GetFullFileName(1,sys,sNN,field));
3e478dba 426 fAcceptance =
427 dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
428 if(!fAcceptance) {
429 ::Error("Init", "no central Acceptance found!") ;
430 return;
431 }
e58000b7 432
433 if(fSecmap && fAcceptance) {
434 fIsInit = kTRUE;
435 ::Info("Init",
436 "Central Manager initialised for sys %d, energy %d, field %d",sys,sNN,field);
6f791cc3 437
e58000b7 438 }
439
6f791cc3 440}
441//
442// EOF
443//