]>
Commit | Line | Data |
---|---|---|
ea3fd2d5 | 1 | // $Id$ |
2 | ||
ea3fd2d5 | 3 | #include "AliAnalysisTaskEMCALPi0PbPb.h" |
fa443410 | 4 | #include <TAxis.h> |
717fe7de | 5 | #include <TChain.h> |
6 | #include <TClonesArray.h> | |
f5d4ab70 | 7 | #include <TFile.h> |
717fe7de | 8 | #include <TH1F.h> |
9 | #include <TH2F.h> | |
10 | #include <TList.h> | |
11 | #include <TLorentzVector.h> | |
f5d4ab70 | 12 | #include <TNtuple.h> |
717fe7de | 13 | #include "AliAODEvent.h" |
14 | #include "AliAODVertex.h" | |
15 | #include "AliAnalysisManager.h" | |
16 | #include "AliAnalysisTaskEMCALClusterizeFast.h" | |
ea3fd2d5 | 17 | #include "AliCentrality.h" |
ea3fd2d5 | 18 | #include "AliEMCALGeoUtils.h" |
717fe7de | 19 | #include "AliESDEvent.h" |
20 | #include "AliESDVertex.h" | |
43cfaa06 | 21 | #include "AliLog.h" |
ea3fd2d5 | 22 | |
23 | ClassImp(AliAnalysisTaskEMCALPi0PbPb) | |
717fe7de | 24 | |
ea3fd2d5 | 25 | //________________________________________________________________________ |
26 | AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name) | |
27 | : AliAnalysisTaskSE(name), | |
d595acbb | 28 | fAsymMax(1), |
717fe7de | 29 | fCentVar("V0M"), |
30 | fCentFrom(0), | |
31 | fCentTo(100), | |
76332037 | 32 | fVtxZMin(-10), |
33 | fVtxZMax(+10), | |
43cfaa06 | 34 | fUseQualFlag(1), |
717fe7de | 35 | fClusName(), |
f5d4ab70 | 36 | fDoNtuple(0), |
d595acbb | 37 | fGeom(0), |
717fe7de | 38 | fOutput(0), |
39 | fEsdEv(0), | |
40 | fAodEv(0), | |
43cfaa06 | 41 | fRecPoints(0), |
717fe7de | 42 | fEsdClusters(0), |
43 | fEsdCells(0), | |
44 | fAodClusters(0), | |
286b47a5 | 45 | fAodCells(0), |
fa443410 | 46 | fPtRanges(0), |
f5d4ab70 | 47 | fNtuple(0), |
fa443410 | 48 | fHCuts(0x0), |
49 | fHVertexZ(0x0), | |
76332037 | 50 | fHVertexZ2(0x0), |
fa443410 | 51 | fHCent(0x0), |
52 | fHCentQual(0x0), | |
d595acbb | 53 | fHColuRow(0x0), |
54 | fHColuRowE(0x0), | |
55 | fHCellMult(0x0), | |
56 | fHCellE(0x0), | |
57 | fHCellH(0x0), | |
fa443410 | 58 | fHClustEtaPhi(0x0), |
59 | fHClustEnergyPt(0x0), | |
60 | fHClustEnergySigma(0x0), | |
d595acbb | 61 | fHClustSigmaSigma(0x0), |
fa443410 | 62 | fHClustNTowEnergyRatio(0x0), |
63 | fHPionEtaPhi(0x0), | |
64 | fHPionMggPt(0x0), | |
65 | fHPionMggAsym(0x0) | |
ea3fd2d5 | 66 | { |
717fe7de | 67 | // Constructor. |
ea3fd2d5 | 68 | |
d595acbb | 69 | if (!name) |
70 | return; | |
71 | SetName(name); | |
ea3fd2d5 | 72 | DefineInput(0, TChain::Class()); |
ea3fd2d5 | 73 | DefineOutput(1, TList::Class()); |
fa443410 | 74 | fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells. " |
29b496d5 | 75 | "AOD:header,vertices,emcalCells"; |
ea3fd2d5 | 76 | } |
ea3fd2d5 | 77 | |
717fe7de | 78 | //________________________________________________________________________ |
79 | AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb() | |
80 | { | |
81 | // Destructor. | |
ea3fd2d5 | 82 | |
717fe7de | 83 | delete fOutput; fOutput = 0; |
fa443410 | 84 | delete fPtRanges; fPtRanges = 0; |
d595acbb | 85 | delete fGeom; fGeom = 0; |
86 | delete [] fHColuRow; | |
87 | delete [] fHColuRowE; | |
88 | delete [] fHCellMult; | |
ea3fd2d5 | 89 | } |
717fe7de | 90 | |
ea3fd2d5 | 91 | //________________________________________________________________________ |
92 | void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects() | |
93 | { | |
717fe7de | 94 | // Create user objects here. |
ea3fd2d5 | 95 | |
717fe7de | 96 | fOutput = new TList(); |
97 | fOutput->SetOwner(); | |
ea3fd2d5 | 98 | |
f5d4ab70 | 99 | if (fDoNtuple) { |
100 | TFile *f = OpenFile(1); | |
101 | if (f) | |
102 | fNtuple = new TNtuple("nt","nt","cent:m:pt:e1:e2:e1m:e2m:n1:n2:db1:db2:disp1:disp2:mm1:mm2:ms1:ms2"); | |
103 | } | |
104 | ||
d595acbb | 105 | // histograms |
fa443410 | 106 | fHCuts = new TH1F("hEventCuts","",4,0.5,4.5); |
107 | fHCuts->GetXaxis()->SetBinLabel(1,"All (PS)"); | |
108 | fHCuts->GetXaxis()->SetBinLabel(2,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo)); | |
109 | fHCuts->GetXaxis()->SetBinLabel(3,"QFlag"); | |
110 | fHCuts->GetXaxis()->SetBinLabel(4,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax)); | |
111 | fOutput->Add(fHCuts); | |
112 | fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25); | |
113 | fHVertexZ->SetXTitle("z [cm]"); | |
114 | fOutput->Add(fHVertexZ); | |
76332037 | 115 | fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25); |
116 | fHVertexZ2->SetXTitle("z [cm]"); | |
117 | fOutput->Add(fHVertexZ2); | |
fa443410 | 118 | fHCent = new TH1F("hCentBeforeCut","",101,-1,100); |
119 | fHCent->SetXTitle(fCentVar.Data()); | |
76332037 | 120 | fOutput->Add(fHCent); |
fa443410 | 121 | fHCentQual = new TH1F("hCentAfterCut","",101,-1,100); |
122 | fHCentQual->SetXTitle(fCentVar.Data()); | |
123 | fOutput->Add(fHCentQual); | |
90d5b88b | 124 | |
d595acbb | 125 | // histograms for cells |
126 | fHColuRow = new TH2F*[4]; | |
127 | fHColuRowE = new TH2F*[4]; | |
128 | fHCellMult = new TH1F*[4]; | |
129 | for (Int_t i = 0; i<4; ++i) { | |
130 | fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",49,0,49,25,0,25); | |
131 | fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i)); | |
132 | fHColuRow[i]->SetXTitle("col (i#eta)"); | |
133 | fHColuRow[i]->SetYTitle("row (i#phi)"); | |
134 | fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",49,0,49,25,0,25); | |
90d5b88b | 135 | fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i)); |
d595acbb | 136 | fHColuRowE[i]->SetXTitle("col (i#eta)"); |
137 | fHColuRowE[i]->SetYTitle("row (i#phi)"); | |
138 | fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000); | |
139 | fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i)); | |
90d5b88b | 140 | fHCellMult[i]->SetXTitle("# of cells"); |
d595acbb | 141 | fOutput->Add(fHColuRow[i]); |
142 | fOutput->Add(fHColuRowE[i]); | |
143 | fOutput->Add(fHCellMult[i]); | |
144 | } | |
145 | fHCellE = new TH1F("hCellE","",100,0.,10.); | |
146 | fHCellE->SetXTitle("E_{cell} [GeV]"); | |
147 | fOutput->Add(fHCellE); | |
148 | fHCellH = new TH1F ("fHCellHighestE","",100,0.,10.); | |
149 | fHCellH->SetXTitle("E^{max}_{cell} [GeV] "); | |
150 | fOutput->Add(fHCellH); | |
90d5b88b | 151 | |
d595acbb | 152 | // histograms for clusters |
fa443410 | 153 | fHClustEtaPhi = new TH2F("hClustEtaPhi","",100,-0.8,0.8,100,1.2,2.2); |
154 | fHClustEtaPhi->SetXTitle("#eta"); | |
155 | fHClustEtaPhi->SetYTitle("#varphi"); | |
156 | fOutput->Add(fHClustEtaPhi); | |
157 | fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50); | |
158 | fHClustEnergyPt->SetXTitle("E [GeV]"); | |
159 | fHClustEnergyPt->SetYTitle("p_{T}^{} [GeV/c]"); | |
160 | fOutput->Add(fHClustEnergyPt); | |
d595acbb | 161 | fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,100,0,30); |
162 | fHClustEnergySigma->SetXTitle("E_{C}^{}*#sigma_{max}^{} [GeV*cm]"); | |
163 | fHClustEnergySigma->SetYTitle("E_{C}^{} [GeV]"); | |
164 | fOutput->Add(fHClustEnergySigma); | |
165 | fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",100,0,100,100,0,20); | |
166 | fHClustSigmaSigma->SetXTitle("#lambda_{0}^{} [cm]"); | |
167 | fHClustSigmaSigma->SetYTitle("#sigma_{max}^{} [cm]"); | |
168 | fOutput->Add(fHClustSigmaSigma); | |
fa443410 | 169 | fHClustNTowEnergyRatio = new TH2F("hClustNTowEnergyRatio","",15,-0.5,14.5,101,-0.05,1.05); |
170 | fHClustNTowEnergyRatio->SetXTitle("N towers"); | |
171 | fHClustNTowEnergyRatio->SetYTitle("Energy ratio"); | |
d595acbb | 172 | fOutput->Add(fHClustNTowEnergyRatio); |
90d5b88b | 173 | |
d595acbb | 174 | // histograms for pion candidates |
fa443410 | 175 | fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100,1.2,2.2); |
176 | fHPionEtaPhi->SetXTitle("#eta^{#gamma#gamma}"); | |
177 | fHPionEtaPhi->SetYTitle("#varphi^{#gamma#gamma}"); | |
178 | fOutput->Add(fHPionEtaPhi); | |
179 | fHPionMggPt = new TH2F("hPionMggPt","",100,0,2,100,0,20.0); | |
180 | fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); | |
181 | fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]"); | |
182 | fOutput->Add(fHPionMggPt); | |
183 | fHPionMggAsym = new TH2F("hPionMggAsym","",100,0,2,100,0,1); | |
184 | fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); | |
185 | fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]"); | |
186 | fOutput->Add(fHPionMggAsym); | |
187 | const Int_t nbins = 19; | |
188 | Double_t xbins[nbins] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,6,7,8,9,10,12.5,15,20,25,50}; | |
189 | fPtRanges = new TAxis(nbins-1,xbins); | |
190 | for (Int_t i = 0; i<=nbins; ++i) { | |
191 | fHPionInvMasses[i] = new TH1F(Form("HPionInvMass%d",i),"",100,0,2); | |
192 | fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); | |
193 | if (i==0) | |
194 | fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0])); | |
195 | else if (i==nbins) | |
196 | fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50")); | |
197 | else | |
198 | fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i])); | |
fa443410 | 199 | fOutput->Add(fHPionInvMasses[i]); |
200 | } | |
d595acbb | 201 | |
717fe7de | 202 | PostData(1, fOutput); |
ea3fd2d5 | 203 | } |
204 | ||
205 | //________________________________________________________________________ | |
206 | void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *) | |
207 | { | |
717fe7de | 208 | // Called for each event. |
209 | ||
210 | if (!InputEvent()) | |
211 | return; | |
212 | ||
43cfaa06 | 213 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); |
214 | fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent()); | |
215 | if (fEsdEv) { | |
216 | am->LoadBranch("AliESDRun."); | |
217 | am->LoadBranch("AliESDHeader."); | |
218 | } else { | |
219 | fAodEv = dynamic_cast<AliAODEvent*>(InputEvent()); | |
220 | am->LoadBranch("header"); | |
221 | } | |
222 | ||
90d5b88b | 223 | if (!fGeom) { // set misalignment matrices (stored in first event) |
d595acbb | 224 | fGeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL"); |
225 | if (fEsdEv) { | |
226 | for (Int_t i=0; i<4; ++i) | |
227 | fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i); | |
228 | } else { | |
229 | for (Int_t i=0; i<4; ++i) | |
230 | fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i); | |
231 | } | |
232 | fGeom->GetEMCGeometry(); | |
233 | } | |
234 | ||
286b47a5 | 235 | Int_t cut = 1; |
fa443410 | 236 | fHCuts->Fill(cut++); |
286b47a5 | 237 | |
717fe7de | 238 | const AliCentrality *centP = InputEvent()->GetCentrality(); |
239 | Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar); | |
fa443410 | 240 | fHCent->Fill(cent); |
717fe7de | 241 | if (cent<fCentFrom||cent>fCentTo) |
286b47a5 | 242 | return; |
43cfaa06 | 243 | |
fa443410 | 244 | fHCuts->Fill(cut++); |
286b47a5 | 245 | |
43cfaa06 | 246 | if (fUseQualFlag) { |
247 | if (centP->GetQuality()>0) | |
248 | return; | |
249 | } | |
286b47a5 | 250 | |
fa443410 | 251 | fHCentQual->Fill(cent); |
252 | fHCuts->Fill(cut++); | |
717fe7de | 253 | |
43cfaa06 | 254 | if (fEsdEv) { |
fa443410 | 255 | am->LoadBranch("PrimaryVertex."); |
256 | am->LoadBranch("SPDVertex."); | |
257 | am->LoadBranch("TPCVertex."); | |
43cfaa06 | 258 | } else { |
259 | fAodEv = dynamic_cast<AliAODEvent*>(InputEvent()); | |
260 | am->LoadBranch("vertices"); | |
261 | } | |
262 | ||
263 | const AliVVertex *vertex = InputEvent()->GetPrimaryVertex(); | |
717fe7de | 264 | if (!vertex) |
265 | return; | |
266 | ||
fa443410 | 267 | fHVertexZ->Fill(vertex->GetZ()); |
286b47a5 | 268 | |
717fe7de | 269 | if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax) |
270 | return; | |
286b47a5 | 271 | |
fa443410 | 272 | fHCuts->Fill(cut++); |
76332037 | 273 | fHVertexZ2->Fill(vertex->GetZ()); |
717fe7de | 274 | |
43cfaa06 | 275 | fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used |
276 | fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached | |
277 | fEsdCells = 0; // will be set if ESD input used | |
278 | fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached | |
279 | // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode | |
280 | fAodCells = 0; // will be set if AOD input used | |
717fe7de | 281 | |
43cfaa06 | 282 | // deal with special output from AliAnalysisTaskEMCALClusterizeFast first |
283 | Bool_t clusattached = 0; | |
284 | Bool_t recalibrated = 0; | |
285 | if (1 && !fClusName.IsNull()) { | |
286 | AliAnalysisTaskEMCALClusterizeFast *cltask = 0; | |
287 | TObjArray *ts = am->GetTasks(); | |
288 | cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName)); | |
289 | if (cltask && cltask->GetClusters()) { | |
d595acbb | 290 | fRecPoints = const_cast<TObjArray*>(cltask->GetClusters()); |
43cfaa06 | 291 | clusattached = cltask->GetAttachClusters(); |
292 | if (cltask->GetCalibData()!=0) | |
293 | recalibrated = kTRUE; | |
294 | } | |
295 | } | |
296 | if (1 && AODEvent() && !fClusName.IsNull()) { | |
717fe7de | 297 | TList *l = AODEvent()->GetList(); |
298 | TClonesArray *clus = 0; | |
299 | if (l) { | |
300 | clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName)); | |
301 | fAodClusters = clus; | |
302 | } | |
ea3fd2d5 | 303 | } |
717fe7de | 304 | |
305 | if (fEsdEv) { // ESD input mode | |
43cfaa06 | 306 | if (1 && (!fRecPoints||clusattached)) { |
307 | if (!clusattached) | |
308 | am->LoadBranch("CaloClusters"); | |
717fe7de | 309 | TList *l = fEsdEv->GetList(); |
310 | TClonesArray *clus = 0; | |
311 | if (l) { | |
312 | clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters")); | |
313 | fEsdClusters = clus; | |
ea3fd2d5 | 314 | } |
315 | } | |
43cfaa06 | 316 | if (1) { |
317 | if (!recalibrated) | |
318 | am->LoadBranch("EMCALCells."); | |
717fe7de | 319 | fEsdCells = fEsdEv->GetEMCALCells(); |
320 | } | |
321 | } else if (fAodEv) { // AOD input mode | |
43cfaa06 | 322 | if (1 && (!fAodClusters || clusattached)) { |
323 | if (!clusattached) | |
324 | am->LoadBranch("caloClusters"); | |
717fe7de | 325 | TList *l = fAodEv->GetList(); |
326 | TClonesArray *clus = 0; | |
327 | if (l) { | |
328 | clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters")); | |
329 | fAodClusters = clus; | |
ea3fd2d5 | 330 | } |
331 | } | |
43cfaa06 | 332 | if (1) { |
333 | if (!recalibrated) | |
334 | am->LoadBranch("emcalCells"); | |
717fe7de | 335 | fAodCells = fAodEv->GetEMCALCells(); |
336 | } | |
337 | } else { | |
338 | AliFatal("Impossible to not have either pointer to ESD or AOD event"); | |
339 | } | |
ea3fd2d5 | 340 | |
43cfaa06 | 341 | if (1) { |
342 | AliDebug(2,Form("fRecPoints set: %p", fRecPoints)); | |
343 | AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters)); | |
344 | AliDebug(2,Form("fEsdCells set: %p", fEsdCells)); | |
345 | AliDebug(2,Form("fAodClusters set: %p", fAodClusters)); | |
346 | AliDebug(2,Form("fAodCells set: %p", fAodCells)); | |
347 | } | |
348 | ||
286b47a5 | 349 | FillCellHists(); |
350 | FillClusHists(); | |
351 | FillPionHists(); | |
ea3fd2d5 | 352 | |
717fe7de | 353 | PostData(1, fOutput); |
ea3fd2d5 | 354 | } |
355 | ||
356 | //________________________________________________________________________ | |
357 | void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *) | |
358 | { | |
717fe7de | 359 | // Terminate called at the end of analysis. |
f5d4ab70 | 360 | |
361 | if (fNtuple) { | |
362 | TFile *f = OpenFile(1); | |
363 | if (f) | |
364 | fNtuple->Write(); | |
365 | } | |
ea3fd2d5 | 366 | } |
ea3fd2d5 | 367 | |
286b47a5 | 368 | //________________________________________________________________________ |
369 | void AliAnalysisTaskEMCALPi0PbPb::FillCellHists() | |
370 | { | |
371 | // Fill histograms related to cell properties. | |
d595acbb | 372 | |
90d5b88b | 373 | AliVCaloCells *cells = fEsdCells; |
374 | if (!cells) | |
375 | cells = fAodCells; | |
376 | ||
377 | if (cells) { | |
378 | Int_t cellModCount[4] = {0,0,0,0}; | |
379 | Double_t cellMaxE = 0; | |
380 | Int_t ncells = cells->GetNumberOfCells(); | |
381 | ||
d595acbb | 382 | for (Int_t i = 0; i<ncells; ++i ) { |
90d5b88b | 383 | Int_t absID = cells->GetCellNumber(i); |
384 | Double_t cellE = cells->GetAmplitude(i); | |
d595acbb | 385 | fHCellE->Fill(cellE); |
90d5b88b | 386 | if (cellE>cellMaxE) |
387 | cellMaxE = cellE; | |
d595acbb | 388 | |
90d5b88b | 389 | Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1; |
390 | Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta); | |
391 | if (!ret) { | |
392 | AliError(Form("Could not get cell index for %d", absID)); | |
393 | continue; | |
394 | } | |
395 | ++cellModCount[iSM]; | |
396 | Int_t iPhi=-1, iEta=-1; | |
397 | fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta); | |
398 | fHColuRow[iSM]->Fill(iEta,iPhi,1); | |
399 | fHColuRowE[iSM]->Fill(iEta,iPhi,cellE); | |
400 | } | |
401 | fHCellH->Fill(cellMaxE); | |
402 | for (Int_t i=0; i<4; ++i) | |
403 | fHCellMult[i]->Fill(cellModCount[i]); | |
404 | } | |
286b47a5 | 405 | } |
406 | ||
407 | //________________________________________________________________________ | |
408 | void AliAnalysisTaskEMCALPi0PbPb::FillClusHists() | |
409 | { | |
90d5b88b | 410 | // Fill histograms related to cluster properties. |
fa443410 | 411 | |
412 | Double_t vertex[3] = {0,0,0}; | |
413 | InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); | |
414 | ||
90d5b88b | 415 | TObjArray *clusters = fEsdClusters; |
416 | if (!clusters) | |
417 | clusters = fAodClusters; | |
418 | ||
419 | if (clusters) { | |
420 | TLorentzVector clusterVec; | |
421 | Int_t nclus = clusters->GetEntries(); | |
fa443410 | 422 | for(Int_t i = 0; i<nclus; ++i) { |
90d5b88b | 423 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); |
fa443410 | 424 | if (!clus) |
425 | continue; | |
426 | if (!clus->IsEMCAL()) | |
427 | continue; | |
428 | clus->GetMomentum(clusterVec,vertex); | |
429 | ||
430 | fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi()); | |
431 | fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt()); | |
d595acbb | 432 | fHClustEnergySigma->Fill(clus->E()*GetSigmaMax(clus),clus->E()); |
433 | fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*GetSigmaMax(clus)); | |
90d5b88b | 434 | fHClustNTowEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)); |
fa443410 | 435 | } |
436 | } | |
286b47a5 | 437 | } |
438 | ||
439 | //________________________________________________________________________ | |
440 | void AliAnalysisTaskEMCALPi0PbPb::FillPionHists() | |
441 | { | |
442 | // Fill histograms related to pions. | |
286b47a5 | 443 | |
fa443410 | 444 | Double_t vertex[3] = {0,0,0}; |
445 | InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); | |
446 | ||
90d5b88b | 447 | TObjArray *clusters = fEsdClusters; |
448 | if (!clusters) | |
449 | clusters = fAodClusters; | |
fa443410 | 450 | |
90d5b88b | 451 | if (clusters) { |
452 | TLorentzVector clusterVec1; | |
453 | TLorentzVector clusterVec2; | |
454 | TLorentzVector pionVec; | |
455 | ||
456 | Int_t nclus = clusters->GetEntries(); | |
fa443410 | 457 | for (Int_t i = 0; i<nclus; ++i) { |
90d5b88b | 458 | AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i)); |
fa443410 | 459 | if (!clus1) |
460 | continue; | |
461 | if (!clus1->IsEMCAL()) | |
462 | continue; | |
463 | clus1->GetMomentum(clusterVec1,vertex); | |
464 | for (Int_t j = i+1; j<nclus; ++j) { | |
90d5b88b | 465 | AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j)); |
fa443410 | 466 | if (!clus2) |
467 | continue; | |
468 | if (!clus2->IsEMCAL()) | |
469 | continue; | |
470 | clus2->GetMomentum(clusterVec2,vertex); | |
471 | pionVec = clusterVec1 + clusterVec2; | |
472 | Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E(); | |
d595acbb | 473 | if (pionZgg < fAsymMax) { |
474 | fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi()); | |
475 | fHPionMggPt->Fill(pionVec.M(),pionVec.Pt()); | |
476 | fHPionMggAsym->Fill(pionVec.M(),pionZgg); | |
477 | Int_t bin = fPtRanges->FindBin(pionVec.Pt()); | |
478 | fHPionInvMasses[bin]->Fill(pionVec.M()); | |
479 | } | |
f5d4ab70 | 480 | |
481 | if (fNtuple) { | |
482 | Double_t mass = pionVec.M(); | |
483 | if (mass>0.08 && mass<0.2) { | |
484 | Float_t vals[17]; | |
485 | vals[0] = InputEvent()->GetCentrality()->GetCentralityPercentileUnchecked(fCentVar); | |
486 | vals[1] = mass; | |
487 | vals[2] = pionVec.Pt(); | |
488 | vals[3] = clusterVec1.E(); | |
489 | vals[4] = clusterVec2.E(); | |
490 | vals[5] = GetMaxCellEnergy(clus1); | |
491 | vals[6] = GetMaxCellEnergy(clus2); | |
492 | vals[7] = clus1->GetNCells(); | |
493 | vals[8] = clus2->GetNCells(); | |
494 | vals[9] = clus1->GetDistanceToBadChannel(); | |
495 | vals[10] = clus2->GetDistanceToBadChannel(); | |
496 | vals[11] = clus1->GetDispersion(); | |
497 | vals[12] = clus2->GetDispersion(); | |
498 | vals[13] = clus1->GetM20(); | |
499 | vals[14] = clus2->GetM20(); | |
500 | vals[15] = clus1->GetM02(); | |
501 | vals[16] = clus2->GetM02(); | |
502 | fNtuple->Fill(vals); | |
503 | } | |
504 | } | |
fa443410 | 505 | } |
506 | } | |
90d5b88b | 507 | } |
fa443410 | 508 | } |
d595acbb | 509 | |
510 | //________________________________________________________________________ | |
90d5b88b | 511 | Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster) |
d595acbb | 512 | { |
90d5b88b | 513 | // Get maximum energy of attached cell. |
514 | ||
515 | Double_t maxe = 0; | |
516 | ||
517 | Int_t ncells = cluster->GetNCells(); | |
518 | if (fEsdCells) { | |
519 | for (Int_t i=0; i<ncells; i++) { | |
520 | Double_t e = fEsdCells->GetCellAmplitude(cluster->GetCellAbsId(i)); | |
521 | if (e>maxe) | |
522 | maxe = e; | |
523 | } | |
524 | } else { | |
525 | for (Int_t i=0; i<ncells; i++) { | |
526 | Double_t e = fAodCells->GetCellAmplitude(cluster->GetCellAbsId(i)); | |
527 | if (e>maxe) | |
528 | maxe = e; | |
529 | } | |
530 | } | |
531 | return maxe; | |
d595acbb | 532 | } |
533 | ||
534 | //________________________________________________________________________ | |
90d5b88b | 535 | Double_t AliAnalysisTaskEMCALPi0PbPb::GetSigmaMax(AliVCluster *cluster) |
d595acbb | 536 | { |
90d5b88b | 537 | // Calculate the (E) weighted variance along the longer (eigen) axis. |
538 | ||
539 | Double_t sigmaMax = 0; // cluster variance along its longer axis | |
540 | Double_t Ec = cluster->E(); // cluster energy | |
541 | Double_t Xc = 0 ; // cluster first moment along X | |
542 | Double_t Yc = 0 ; // cluster first moment along Y | |
543 | Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2) | |
544 | Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2) | |
545 | Double_t Syy = 0 ; // cluster covariance^2 | |
546 | ||
f5d4ab70 | 547 | Double_t testenergy = 0; |
90d5b88b | 548 | AliVCaloCells *cells = fEsdCells; |
549 | if (!cells) | |
550 | cells = fAodCells; | |
551 | ||
552 | if (cells) { | |
d595acbb | 553 | TVector3 pos; |
90d5b88b | 554 | Int_t ncells = cluster->GetNCells(); |
f5d4ab70 | 555 | if (ncells==1) |
556 | return 0; | |
90d5b88b | 557 | for(Int_t j=0; j<ncells; ++j) { |
558 | fGeom->GetGlobal(cluster->GetCellAbsId(j),pos); | |
559 | Int_t id = cluster->GetCellAbsId(j); | |
f5d4ab70 | 560 | Double_t cellen = cells->GetCellAmplitude(id); |
561 | testenergy += cellen; | |
562 | Xc += cellen*pos.X(); | |
563 | Yc += cellen*pos.Y(); | |
564 | Sxx += cellen*pos.X()*pos.X(); | |
565 | Syy += cellen*pos.Y()*pos.Y(); | |
566 | Sxy += cellen*pos.X()*pos.Y(); | |
90d5b88b | 567 | } |
568 | Xc /= Ec; | |
569 | Yc /= Ec; | |
570 | Sxx /= Ec; | |
571 | Syy /= Ec; | |
572 | Sxy /= Ec; | |
573 | Sxx -= Xc*Xc; | |
574 | Syy -= Yc*Yc; | |
575 | Sxy -= Xc*Yc; | |
576 | sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0; | |
577 | sigmaMax = TMath::Sqrt(sigmaMax); | |
f5d4ab70 | 578 | |
579 | printf("%f %f\n",testenergy,cluster->E()); | |
580 | cout << testenergy << " " << cluster->E() << endl; | |
d595acbb | 581 | } |
d595acbb | 582 | return sigmaMax; |
583 | } | |
90d5b88b | 584 |