1 //====================================================================
3 // Base class for classes that calculate the multiplicity in the
4 // central region event-by-event
10 // - AliAODCentralMult
15 #include "AliCentralMultiplicityTask.h"
16 #include "AliForwardCorrectionManager.h"
17 #include "AliForwardUtil.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"
25 #include "AliFMDEventInspector.h"
32 //====================================================================
33 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
34 : AliAnalysisTaskSE(name),
41 fFirstEventSeen(false)
46 DefineOutput(1, TList::Class());
48 "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
49 "SPDVertex.,PrimaryVertex.";
51 //____________________________________________________________________
52 AliCentralMultiplicityTask::AliCentralMultiplicityTask()
53 : AliAnalysisTaskSE(),
60 fFirstEventSeen(false)
66 //____________________________________________________________________
67 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
68 : AliAnalysisTaskSE(o),
71 fAODCentral(o.fAODCentral),
73 fUseSecondary(o.fUseSecondary),
74 fUseAcceptance(o.fUseAcceptance),
75 fFirstEventSeen(o.fFirstEventSeen)
81 //____________________________________________________________________
82 AliCentralMultiplicityTask&
83 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
86 // Assignment operator
90 fAODCentral = o.fAODCentral;
91 fManager = o.fManager;
92 fUseSecondary = o.fUseSecondary;
93 fUseAcceptance = o.fUseAcceptance;
94 fFirstEventSeen = o.fFirstEventSeen;
97 //____________________________________________________________________
98 void AliCentralMultiplicityTask::UserCreateOutputObjects()
101 // Create output objects
105 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
107 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
108 if (!ah) AliFatal("No AOD output handler set in analysis manager");
111 TObject* obj = &fAODCentral;
112 ah->AddBranch("AliAODCentralMult", &obj);
119 //____________________________________________________________________
120 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
123 // Process each event
129 AliESDInputHandler* eventHandler =
130 dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
131 ->GetInputEventHandler());
133 AliWarning("No inputhandler found for this event!");
136 AliESDEvent* esd = eventHandler->GetEvent();
138 if(!GetManager().IsInit() && !fFirstEventSeen) {
139 AliFMDEventInspector inspector;
140 inspector.ReadRunDetails(esd);
141 GetManager().Init(inspector.GetCollisionSystem(),
142 inspector.GetEnergy(),
143 inspector.GetField());
145 AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
146 fFirstEventSeen = kTRUE;
149 //Selecting only events with |valid vertex| < 10 cm
150 const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
152 if(!(vertex->GetStatus())) return;
153 if(vertex->GetNContributors() <= 0) return ;
154 if(vertex->GetZRes() > 0.1 ) return;
155 Double_t vertexXYZ[3]={0,0,0};
156 vertex->GetXYZ(vertexXYZ);
157 if(TMath::Abs(vertexXYZ[2]) > 10) return;
160 Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
161 //HHD: The vtxbins are 1-10, not 0-9
162 Int_t vtxbin = Int_t(vertexBinDouble + 1) ;
164 // Make sure AOD is filled
165 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
167 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
169 AliFatal("No AOD output handler set in analysis manager");
171 ah->SetFillAOD(kTRUE);
174 fAODCentral.Clear("");
175 TH2D *aodHist = &(fAODCentral.GetHistogram());
177 const AliMultiplicity* spdmult = esd->GetMultiplicity();
178 //Filling clusters in layer 1 used for tracklets...
179 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
180 aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
182 //...and then the unused ones in layer 1
183 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
184 aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
185 spdmult->GetPhiSingle(j));
189 TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
190 TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
192 if (!hSecMap) AliFatal("No secondary map!");
193 if (!hAcceptance) AliFatal("No acceptance!");
195 if (fUseSecondary && hSecMap) aodHist->Divide(hSecMap);
197 for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
198 Float_t accCor = hAcceptance->GetBinContent(nx);
199 Float_t accErr = hAcceptance->GetBinError(nx);
201 Bool_t etabinSeen = kFALSE;
202 for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
203 // Get currrent value
204 Float_t aodValue = aodHist->GetBinContent(nx,ny);
205 Float_t aodErr = aodHist->GetBinError(nx,ny);
209 if(hSecMap) secCor = hSecMap->GetBinContent(nx,ny);
210 if (secCor > 0.5) etabinSeen = kTRUE;
211 if (aodValue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
213 if (!fUseAcceptance) continue;
215 // Acceptance correction
216 if (accCor < 0.000001) accCor = 1;
217 Float_t aodNew = aodValue / accCor ;
218 Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
219 TMath::Power(accErr/accCor,2) );
220 aodHist->SetBinContent(nx,ny, aodNew);
222 aodHist->SetBinError(nx,ny,error);
223 aodHist->SetBinError(nx,ny,aodErr);
226 //Filling underflow bin if we eta bin is in range
227 if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
232 //____________________________________________________________________
233 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
242 //____________________________________________________________________
244 AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
253 //====================================================================
254 AliCentralMultiplicityTask::Manager::Manager() :
255 fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
256 fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
259 fAcceptanceName("centralacceptance"),
260 fSecMapName("centralsecmap"),
267 //____________________________________________________________________
268 AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
269 :fAcceptancePath(o.fAcceptancePath),
270 fSecMapPath(o.fSecMapPath),
271 fAcceptance(o.fAcceptance),
273 fAcceptanceName(o.fAcceptanceName),
274 fSecMapName(o.fSecMapName),
281 //____________________________________________________________________
282 AliCentralMultiplicityTask::Manager&
283 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
286 // Assignment operator
288 fAcceptancePath = o.fAcceptancePath;
289 fSecMapPath = o.fSecMapPath;
290 fAcceptance = o.fAcceptance;
292 fAcceptanceName = o.fAcceptanceName;
293 fSecMapName = o.fSecMapName;
298 //____________________________________________________________________
300 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
306 // Get full path name to object file
310 // sys Collision system
311 // sNN Center of mass energy
312 // field Magnetic field
318 what == 0 ? GetSecMapPath() : GetAcceptancePath(),
319 GetFileName(what, sys, sNN, field));
322 //____________________________________________________________________
324 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
330 // Get the full path name
334 // sys Collision system
335 // sNN Center of mass energy
336 // field Magnetic field
341 // Must be static - otherwise the data may disappear on return from
342 // this member function
343 static TString fname = "";
347 case 0: fname.Append(fSecMapName.Data()); break;
348 case 1: fname.Append(fAcceptanceName.Data()); break;
350 ::Error("GetFileName",
351 "Invalid indentifier %d for central object, must be 0 or 1!", what);
354 fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
355 AliForwardUtil::CollisionSystemString(sys),
356 sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
361 //____________________________________________________________________
363 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
366 // Get the secondary map
375 ::Warning("GetSecMapCorrection","No secondary map defined");
378 return fSecmap->GetCorrection(vtxbin);
380 //____________________________________________________________________
382 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
386 // Get the acceptance correction
395 ::Warning("GetAcceptanceCorrection","No acceptance map defined");
398 return fAcceptance->GetCorrection(vtxbin);
401 //____________________________________________________________________
403 AliCentralMultiplicityTask::Manager::Init(UShort_t sys,
411 // sys Collision system (1: pp, 2: PbPb)
412 // sNN Center of mass energy per nucleon pair [GeV]
413 // field Magnetic field [kG]
415 if(fIsInit) ::Warning("Init","Already initialised - overriding...");
417 TFile fsec(GetFullFileName(0,sys,sNN,field));
419 dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
421 ::Error("Init", "no central Secondary Map found!") ;
424 TFile facc(GetFullFileName(1,sys,sNN,field));
426 dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
428 ::Error("Init", "no central Acceptance found!") ;
432 if(fSecmap && fAcceptance) {
435 "Central Manager initialised for sys %d, energy %d, field %d",sys,sNN,field);