- make the filling of V0 branch optional
[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),
40 firstEventSeen(false)
6f791cc3 41{
42
43 DefineOutput(1, TList::Class());
44}
45//____________________________________________________________________
9c825779 46AliCentralMultiplicityTask::AliCentralMultiplicityTask()
47 : AliAnalysisTaskSE(),
48 fData(0),
49 fList(0),
50 fAODCentral(),
3b2bfb07 51 fManager(),
e58000b7 52 fUseSecondary(true),
53 firstEventSeen(false)
9c825779 54{
55}
56//____________________________________________________________________
57AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
58 : AliAnalysisTaskSE(o),
59 fData(o.fData),
60 fList(o.fList),
61 fAODCentral(o.fAODCentral),
3b2bfb07 62 fManager(o.fManager),
e58000b7 63 fUseSecondary(o.fUseSecondary),
64 firstEventSeen(o.firstEventSeen)
9c825779 65{
66}
67//____________________________________________________________________
68AliCentralMultiplicityTask&
69AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
70{
e58000b7 71 fData = o.fData;
72 fList = o.fList;
73 fAODCentral = o.fAODCentral;
74 fManager = o.fManager;
75 fUseSecondary = o.fUseSecondary;
76 firstEventSeen = o.firstEventSeen;
9c825779 77 return *this;
78}
79//____________________________________________________________________
3e478dba 80void AliCentralMultiplicityTask::UserCreateOutputObjects()
81{
6f791cc3 82
83 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
84 AliAODHandler* ah =
85 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
86 if (!ah) AliFatal("No AOD output handler set in analysis manager");
87
88
89 TObject* obj = &fAODCentral;
90 ah->AddBranch("AliAODCentralMult", &obj);
91
92 fList = new TList();
9d05ffeb 93 fList->SetOwner();
6f791cc3 94 PostData(1,fList);
95
96}
97//____________________________________________________________________
3e478dba 98void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
99{
6f791cc3 100
101 AliESDInputHandler* eventHandler =
102 dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
103 ->GetInputEventHandler());
104 if (!eventHandler) {
105 AliWarning("No inputhandler found for this event!");
106 return;}
107
108 AliESDEvent* esd = eventHandler->GetEvent();
109
e58000b7 110 if(!GetManager().IsInit() && !firstEventSeen) {
111 AliFMDEventInspector inspector;
112 inspector.ReadRunDetails(esd);
113 GetManager().Init(inspector.GetCollisionSystem(),
114 inspector.GetEnergy(),
115 inspector.GetField());
116
117 //std::cout<<inspector.GetCollisionSystem()<<" "<<inspector.GetEnergy()<<" "<<inspector.GetField()<<std::endl;
118 AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
119 firstEventSeen = kTRUE;
120 }
121
6f791cc3 122 //Selecting only events with |valid vertex| < 10 cm
123 const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
124 if (!vertex) return;
125 if(!(vertex->GetStatus())) return;
126 if(vertex->GetNContributors() <= 0) return ;
127 if(vertex->GetZRes() > 0.1 ) return;
128 Double_t vertexXYZ[3]={0,0,0};
129 vertex->GetXYZ(vertexXYZ);
130 if(TMath::Abs(vertexXYZ[2]) > 10) return;
131
3e478dba 132 Double_t delta = 2 ;
6f791cc3 133 Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
3e478dba 134 //HHD: The vtxbins are 1-10, not 0-9
135 Int_t vtxbin = Int_t(vertexBinDouble + 1) ;
6f791cc3 136
137 // Make sure AOD is filled
138 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
139 AliAODHandler* ah =
140 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
141 if (!ah)
142 AliFatal("No AOD output handler set in analysis manager");
143
144 ah->SetFillAOD(kTRUE);
145
146 //Doing analysis
147 fAODCentral.Clear("");
148 TH2D *aodHist = &(fAODCentral.GetHistogram());
149
150 const AliMultiplicity* spdmult = esd->GetMultiplicity();
151 //Filling clusters in layer 1 used for tracklets...
3e478dba 152 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
6f791cc3 153 aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
3e478dba 154
6f791cc3 155 //...and then the unused ones in layer 1
3e478dba 156 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
6f791cc3 157 aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
3e478dba 158 spdmult->GetPhiSingle(j));
6f791cc3 159
6f791cc3 160
3e478dba 161 // Corrections
3b2bfb07 162 TH2D* hSecMap = 0;
6f791cc3 163 TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
3b2bfb07 164 if (fUseSecondary) hSecMap = fManager.GetSecMapCorrection(vtxbin);
165 if (fUseSecondary && !hSecMap) AliFatal("No secondary map!");
166 if (!hAcceptance) AliFatal("No acceptance!");
dbe8d9ed 167
3b2bfb07 168 if (hSecMap) aodHist->Divide(hSecMap);
6f791cc3 169
170 for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
3b2bfb07 171 Float_t accCor = hAcceptance->GetBinContent(nx);
6f791cc3 172
173 Bool_t etabinSeen = kFALSE;
174 for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
3b2bfb07 175 Float_t aodValue = aodHist->GetBinContent(nx,ny);
176 Float_t secCor = hSecMap->GetBinContent(nx,ny);
177 if (secCor > 0.5) etabinSeen = kTRUE;
178 if (aodValue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
179 if (accCor < 0.000001) accCor = 1;
180 Float_t aodNew = aodValue / accCor ;
181 aodHist->SetBinContent(nx,ny, aodNew);
3e478dba 182 Float_t aodErr = aodHist->GetBinError(nx,ny);
183 Float_t accErr = hAcceptance->GetBinError(nx);
3b2bfb07 184 Float_t error = aodNew *TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
185 TMath::Power(accErr/accCor,2) );
6f791cc3 186 aodHist->SetBinError(nx,ny,error);
187
188 }
189 //Filling underflow bin if we eta bin is in range
190 if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
6f791cc3 191 }
192
193 PostData(1,fList);
6f791cc3 194}
195//____________________________________________________________________
3e478dba 196void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
197{
6f791cc3 198}
199//____________________________________________________________________
200void
201AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
202{
6f791cc3 203}
3e478dba 204//====================================================================
6f791cc3 205AliCentralMultiplicityTask::Manager::Manager() :
206 fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
207 fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
208 fAcceptance(),
209 fSecmap(),
210 fAcceptanceName("centralacceptance"),
e58000b7 211 fSecMapName("centralsecmap"),
212 fIsInit(kFALSE)
6f791cc3 213{
214
215
216}
217//____________________________________________________________________
3e478dba 218AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
219 :fAcceptancePath(o.fAcceptancePath),
220 fSecMapPath(o.fSecMapPath),
221 fAcceptance(o.fAcceptance),
222 fSecmap(o.fSecmap),
223 fAcceptanceName(o.fAcceptanceName),
e58000b7 224 fSecMapName(o.fSecMapName),
225 fIsInit(o.fIsInit)
3e478dba 226{}
227//____________________________________________________________________
228AliCentralMultiplicityTask::Manager&
229AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
230{
231 fAcceptancePath = o.fAcceptancePath;
232 fSecMapPath = o.fSecMapPath;
233 fAcceptance = o.fAcceptance;
234 fSecmap = o.fSecmap;
235 fAcceptanceName = o.fAcceptanceName;
236 fSecMapName = o.fSecMapName;
e58000b7 237 fIsInit = o.fIsInit;
3e478dba 238 return *this;
239}
240
241//____________________________________________________________________
242const char*
243AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
244 UShort_t sys,
245 UShort_t sNN,
246 Short_t field) const
247{
248 return Form("%s/%s",
249 what == 0 ? GetSecMapPath() : GetAcceptancePath(),
250 GetFileName(what, sys, sNN, field));
251}
252
253//____________________________________________________________________
254const char*
255AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
256 UShort_t sys,
257 UShort_t sNN,
258 Short_t field) const
259{
260 // Must be static - otherwise the data may disappear on return from
261 // this member function
262 static TString fname = "";
6f791cc3 263
dbe8d9ed 264 fname = "";
6f791cc3 265 switch(what) {
3e478dba 266 case 0: fname.Append(fSecMapName.Data()); break;
267 case 1: fname.Append(fAcceptanceName.Data()); break;
6f791cc3 268 default:
3e478dba 269 ::Error("GetFileName",
270 "Invalid indentifier %d for central object, must be 0 or 1!", what);
6f791cc3 271 break;
272 }
273 fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
274 AliForwardUtil::CollisionSystemString(sys),
275 sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
276
277 return fname.Data();
6f791cc3 278}
3e478dba 279
6f791cc3 280//____________________________________________________________________
3e478dba 281TH2D*
282AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
283{
284 if (!fSecmap) {
285 ::Warning("GetSecMapCorrection","No secondary map defined");
286 return 0;
287 }
288 return fSecmap->GetCorrection(vtxbin);
289}
290//____________________________________________________________________
291TH1D*
292AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
293 const
294{
295 if (!fAcceptance) {
296 ::Warning("GetAcceptanceCorrection","No acceptance map defined");
297 return 0;
298 }
299 return fAcceptance->GetCorrection(vtxbin);
300}
301
302//____________________________________________________________________
303void
304AliCentralMultiplicityTask::Manager::Init(UShort_t sys,
305 UShort_t sNN,
306 Short_t field)
307{
e58000b7 308 if(fIsInit) ::Warning("Init","Already initialised - overriding...");
309
6f791cc3 310 TFile fsec(GetFullFileName(0,sys,sNN,field));
3e478dba 311 fSecmap =
312 dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
6f791cc3 313 if(!fSecmap) {
3e478dba 314 ::Error("Init", "no central Secondary Map found!") ;
6f791cc3 315 return;
316 }
317 TFile facc(GetFullFileName(1,sys,sNN,field));
3e478dba 318 fAcceptance =
319 dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
320 if(!fAcceptance) {
321 ::Error("Init", "no central Acceptance found!") ;
322 return;
323 }
e58000b7 324
325 if(fSecmap && fAcceptance) {
326 fIsInit = kTRUE;
327 ::Info("Init",
328 "Central Manager initialised for sys %d, energy %d, field %d",sys,sNN,field);
6f791cc3 329
e58000b7 330 }
331
6f791cc3 332}
333//
334// EOF
335//