Changed ranges (Chiara)
[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"
25#include <TROOT.h>
26#include <TFile.h>
27#include <iostream>
28#include <iomanip>
29
30//====================================================================
31AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
32 : AliAnalysisTaskSE(name),
33 fData(0),
34 fList(0),
35 fAODCentral(kFALSE),
36 fManager()
37{
38
39 DefineOutput(1, TList::Class());
40}
41//____________________________________________________________________
42void AliCentralMultiplicityTask::UserCreateOutputObjects() {
43
44 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
45 AliAODHandler* ah =
46 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
47 if (!ah) AliFatal("No AOD output handler set in analysis manager");
48
49
50 TObject* obj = &fAODCentral;
51 ah->AddBranch("AliAODCentralMult", &obj);
52
53 fList = new TList();
54 PostData(1,fList);
55
56}
57//____________________________________________________________________
58void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/) {
59
60 AliESDInputHandler* eventHandler =
61 dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
62 ->GetInputEventHandler());
63 if (!eventHandler) {
64 AliWarning("No inputhandler found for this event!");
65 return;}
66
67 AliESDEvent* esd = eventHandler->GetEvent();
68
69 //Selecting only events with |valid vertex| < 10 cm
70 const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
71 if (!vertex) return;
72 if(!(vertex->GetStatus())) return;
73 if(vertex->GetNContributors() <= 0) return ;
74 if(vertex->GetZRes() > 0.1 ) return;
75 Double_t vertexXYZ[3]={0,0,0};
76 vertex->GetXYZ(vertexXYZ);
77 if(TMath::Abs(vertexXYZ[2]) > 10) return;
78
79 Double_t delta = 2 ;
80 Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
81 Int_t vtxbin = (Int_t)vertexBinDouble +1 ; //HHD: The vtxbins are 1-10, not 0-9
82
83 // Make sure AOD is filled
84 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
85 AliAODHandler* ah =
86 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
87 if (!ah)
88 AliFatal("No AOD output handler set in analysis manager");
89
90 ah->SetFillAOD(kTRUE);
91
92 //Doing analysis
93 fAODCentral.Clear("");
94 TH2D *aodHist = &(fAODCentral.GetHistogram());
95
96 const AliMultiplicity* spdmult = esd->GetMultiplicity();
97 //Filling clusters in layer 1 used for tracklets...
98 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
99 aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
100 }
101 //...and then the unused ones in layer 1
102 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
103 aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
104 spdmult->GetPhiSingle(j));
105 }
106
107 // Corrections
108
109 TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
110 TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
111
112 aodHist->Divide(hSecMap);
113
114 for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
115 Float_t acccor = hAcceptance->GetBinContent(nx);
116
117 Bool_t etabinSeen = kFALSE;
118 for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
119 Float_t aodvalue = aodHist->GetBinContent(nx,ny);
120 Float_t seccor = hSecMap->GetBinContent(nx,ny);
121 if(seccor > 0.5) etabinSeen = kTRUE;
122 if(aodvalue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
123
124 Float_t aodnew = aodvalue / acccor ;
125 aodHist->SetBinContent(nx,ny, aodnew);
126 Float_t error = aodnew*TMath::Sqrt(TMath::Power(aodHist->GetBinError(nx,ny)/aodvalue,2) +
127 TMath::Power(hAcceptance->GetBinError(nx)/acccor,2) );
128 aodHist->SetBinError(nx,ny,error);
129
130 }
131 //Filling underflow bin if we eta bin is in range
132 if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
133
134 }
135
136 PostData(1,fList);
137
138}
139//____________________________________________________________________
140void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/) {
141
142
143
144}
145//____________________________________________________________________
146void
147AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
148{
149
150}
151//____________________________________________________________________
152AliCentralMultiplicityTask::Manager::Manager() :
153 fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
154 fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
155 fAcceptance(),
156 fSecmap(),
157 fAcceptanceName("centralacceptance"),
158 fSecMapName("centralsecmap")
159{
160
161
162}
163//____________________________________________________________________
164const char* AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
165 UShort_t sys,
166 UShort_t sNN,
167 Short_t field) {
168
169 TString fname = "";
170
171 switch(what) {
172 case 0:
173 fname.Append(fSecMapName.Data());
174 break;
175 case 1:
176 fname.Append(fAcceptanceName.Data());
177 break;
178 default:
179 printf("Invalid indentifier for central object, must be 0 or 1!");
180 break;
181 }
182 fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
183 AliForwardUtil::CollisionSystemString(sys),
184 sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
185
186 return fname.Data();
187
188}
189//____________________________________________________________________
190void AliCentralMultiplicityTask::Manager::Init( UShort_t sys,
191 UShort_t sNN,
192 Short_t field) {
193
194 TFile fsec(GetFullFileName(0,sys,sNN,field));
195 fSecmap = dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
196 if(!fSecmap) {
197 printf("no central Secondary Map found!") ;
198 return;
199 }
200 TFile facc(GetFullFileName(1,sys,sNN,field));
201 fAcceptance = dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
202if(!fAcceptance) {
203 printf("no central Acceptance found!") ;
204 return;
205 }
206
207}
208//
209// EOF
210//