]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/UserTasks/EmcalTasks/AliAnalysisTaskEMCALPi0PbPb.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes for #86595 commit in STEER and port to Release AliESDTZERO with fixed coverity
[u/mrichter/AliRoot.git] / PWG4 / UserTasks / EmcalTasks / AliAnalysisTaskEMCALPi0PbPb.cxx
CommitLineData
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>
296ea9b4 8#include <TGeoGlobalMagField.h>
717fe7de 9#include <TH1F.h>
10#include <TH2F.h>
11#include <TList.h>
12#include <TLorentzVector.h>
f5d4ab70 13#include <TNtuple.h>
296ea9b4 14#include <TProfile.h>
b3ee6797 15#include <TString.h>
296ea9b4 16#include <TVector2.h>
717fe7de 17#include "AliAODEvent.h"
56fd6cb2 18#include "AliAODMCParticle.h"
717fe7de 19#include "AliAODVertex.h"
20#include "AliAnalysisManager.h"
21#include "AliAnalysisTaskEMCALClusterizeFast.h"
296ea9b4 22#include "AliCDBManager.h"
ea3fd2d5 23#include "AliCentrality.h"
c2fe5f0e 24#include "AliEMCALDigit.h"
2ef5608f 25#include "AliEMCALGeometry.h"
0fbe8d4f 26#include "AliEMCALRecPoint.h"
296ea9b4 27#include "AliEMCALRecoUtils.h"
3a952328 28#include "AliESDCaloTrigger.h"
717fe7de 29#include "AliESDEvent.h"
5fe1ca23 30#include "AliESDUtils.h"
717fe7de 31#include "AliESDVertex.h"
296ea9b4 32#include "AliESDtrack.h"
0ec74551 33#include "AliESDtrackCuts.h"
b6c599fe 34#include "AliEventplane.h"
296ea9b4 35#include "AliGeomManager.h"
b3ee6797 36#include "AliInputEventHandler.h"
43cfaa06 37#include "AliLog.h"
38727e64 38#include "AliMCEvent.h"
39#include "AliMCParticle.h"
296ea9b4 40#include "AliMagF.h"
5fe1ca23 41#include "AliMultiplicity.h"
38727e64 42#include "AliStack.h"
0ec74551 43#include "AliTrackerBase.h"
c0d2e1f2 44#include "AliTriggerAnalysis.h"
ea3fd2d5 45
46ClassImp(AliAnalysisTaskEMCALPi0PbPb)
717fe7de 47
4ea96211 48//________________________________________________________________________
49AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb()
50 : AliAnalysisTaskSE(),
51 fCentVar("V0M"),
52 fCentFrom(0),
53 fCentTo(100),
54 fVtxZMin(-10),
55 fVtxZMax(+10),
56 fUseQualFlag(1),
57 fClusName(),
58 fDoNtuple(0),
59 fDoAfterburner(0),
60 fAsymMax(1),
61 fNminCells(1),
62 fMinE(0.100),
63 fMinErat(0),
64 fMinEcc(0),
65 fGeoName("EMCAL_FIRSTYEARV1"),
66 fMinNClusPerTr(50),
67 fIsoDist(0.2),
68 fTrClassNames(""),
69 fTrCuts(0),
70 fPrimTrCuts(0),
71 fDoTrMatGeom(0),
72 fTrainMode(0),
73 fMarkCells(),
74 fMinL0Time(-1),
75 fMaxL0Time(1024),
76 fMcMode(0),
77 fEmbedMode(0),
78 fGeom(0),
79 fReco(0),
80 fDoPSel(kTRUE),
81 fIsGeoMatsSet(0),
82 fNEvs(0),
83 fOutput(0),
84 fTrClassNamesArr(0),
85 fEsdEv(0),
86 fAodEv(0),
87 fRecPoints(0),
88 fDigits(0),
89 fEsdClusters(0),
90 fEsdCells(0),
91 fAodClusters(0),
92 fAodCells(0),
93 fPtRanges(0),
94 fSelTracks(0),
95 fSelPrimTracks(0),
96 fNAmpInTrigger(0),
97 fAmpInTrigger(0),
98 fNtuple(0),
99 fHeader(0),
100 fPrimVert(0),
101 fSpdVert(0),
102 fTpcVert(0),
103 fClusters(0),
104 fTriggers(0),
105 fMcParts(0),
106 fHCuts(0x0),
107 fHVertexZ(0x0),
108 fHVertexZ2(0x0),
109 fHCent(0x0),
110 fHCentQual(0x0),
111 fHTclsBeforeCuts(0x0),
112 fHTclsAfterCuts(0x0),
113 fHColuRow(0x0),
114 fHColuRowE(0x0),
115 fHCellMult(0x0),
116 fHCellE(0x0),
117 fHCellH(0x0),
118 fHCellM(0x0),
119 fHCellM2(0x0),
120 fHCellFreqNoCut(0x0),
121 fHCellFreqCut100M(0x0),
122 fHCellFreqCut300M(0x0),
123 fHCellFreqE(0x0),
124 fHCellCheckE(0x0),
125 fHClustEccentricity(0),
126 fHClustEtaPhi(0x0),
127 fHClustEnergyPt(0x0),
128 fHClustEnergySigma(0x0),
129 fHClustSigmaSigma(0x0),
130 fHClustNCellEnergyRatio(0x0),
131 fHClustEnergyNCell(0x0),
132 fHPrimTrackPt(0x0),
133 fHPrimTrackEta(0x0),
134 fHPrimTrackPhi(0x0),
135 fHMatchDr(0x0),
136 fHMatchDz(0x0),
137 fHMatchEp(0x0),
138 fHPionEtaPhi(0x0),
139 fHPionMggPt(0x0),
140 fHPionMggAsym(0x0),
141 fHPionMggDgg(0x0)
142{
143 // Constructor.
144}
145
ea3fd2d5 146//________________________________________________________________________
147AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
148 : AliAnalysisTaskSE(name),
717fe7de 149 fCentVar("V0M"),
150 fCentFrom(0),
151 fCentTo(100),
76332037 152 fVtxZMin(-10),
153 fVtxZMax(+10),
43cfaa06 154 fUseQualFlag(1),
717fe7de 155 fClusName(),
f5d4ab70 156 fDoNtuple(0),
a49742b5 157 fDoAfterburner(0),
f224d35b 158 fAsymMax(1),
a49742b5 159 fNminCells(1),
3c661da5 160 fMinE(0.100),
f224d35b 161 fMinErat(0),
162 fMinEcc(0),
6bf90832 163 fGeoName("EMCAL_FIRSTYEARV1"),
b6c599fe 164 fMinNClusPerTr(50),
296ea9b4 165 fIsoDist(0.2),
b3ee6797 166 fTrClassNames(""),
0ec74551 167 fTrCuts(0),
b6c599fe 168 fPrimTrCuts(0),
169 fDoTrMatGeom(0),
3a952328 170 fTrainMode(0),
171 fMarkCells(),
172 fMinL0Time(-1),
173 fMaxL0Time(1024),
38727e64 174 fMcMode(0),
cfd7d5b2 175 fEmbedMode(0),
d595acbb 176 fGeom(0),
296ea9b4 177 fReco(0),
807016ea 178 fDoPSel(kTRUE),
38727e64 179 fIsGeoMatsSet(0),
180 fNEvs(0),
717fe7de 181 fOutput(0),
b3ee6797 182 fTrClassNamesArr(0),
717fe7de 183 fEsdEv(0),
184 fAodEv(0),
43cfaa06 185 fRecPoints(0),
c2fe5f0e 186 fDigits(0),
717fe7de 187 fEsdClusters(0),
188 fEsdCells(0),
189 fAodClusters(0),
286b47a5 190 fAodCells(0),
fa443410 191 fPtRanges(0),
296ea9b4 192 fSelTracks(0),
b6c599fe 193 fSelPrimTracks(0),
3a952328 194 fNAmpInTrigger(0),
195 fAmpInTrigger(0),
788ca675 196 fNtuple(0),
197 fHeader(0),
198 fPrimVert(0),
199 fSpdVert(0),
200 fTpcVert(0),
201 fClusters(0),
3a952328 202 fTriggers(0),
807016ea 203 fMcParts(0),
fa443410 204 fHCuts(0x0),
205 fHVertexZ(0x0),
76332037 206 fHVertexZ2(0x0),
fa443410 207 fHCent(0x0),
208 fHCentQual(0x0),
b3ee6797 209 fHTclsBeforeCuts(0x0),
210 fHTclsAfterCuts(0x0),
d595acbb 211 fHColuRow(0x0),
212 fHColuRowE(0x0),
213 fHCellMult(0x0),
214 fHCellE(0x0),
215 fHCellH(0x0),
6eb6260e 216 fHCellM(0x0),
a49742b5 217 fHCellM2(0x0),
296ea9b4 218 fHCellFreqNoCut(0x0),
2e4d8148 219 fHCellFreqCut100M(0x0),
220 fHCellFreqCut300M(0x0),
221 fHCellFreqE(0x0),
296ea9b4 222 fHCellCheckE(0x0),
6eb6260e 223 fHClustEccentricity(0),
fa443410 224 fHClustEtaPhi(0x0),
225 fHClustEnergyPt(0x0),
226 fHClustEnergySigma(0x0),
d595acbb 227 fHClustSigmaSigma(0x0),
6eb6260e 228 fHClustNCellEnergyRatio(0x0),
f5e0f1e2 229 fHClustEnergyNCell(0x0),
230 fHPrimTrackPt(0x0),
231 fHPrimTrackEta(0x0),
232 fHPrimTrackPhi(0x0),
b6c599fe 233 fHMatchDr(0x0),
234 fHMatchDz(0x0),
235 fHMatchEp(0x0),
fa443410 236 fHPionEtaPhi(0x0),
237 fHPionMggPt(0x0),
6eb6260e 238 fHPionMggAsym(0x0),
239 fHPionMggDgg(0x0)
ea3fd2d5 240{
717fe7de 241 // Constructor.
ea3fd2d5 242
ea3fd2d5 243 DefineOutput(1, TList::Class());
c0d2e1f2 244 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks,EMCALTrigger. "
296ea9b4 245 "AOD:header,vertices,emcalCells,tracks";
ea3fd2d5 246}
ea3fd2d5 247
717fe7de 248//________________________________________________________________________
249AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
250{
251 // Destructor.
ea3fd2d5 252
b3ee6797 253 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
254 delete fOutput; fOutput = 0;
255 }
fa443410 256 delete fPtRanges; fPtRanges = 0;
a2de5ca1 257 fGeom = 0; // do not delete geometry when using instance
296ea9b4 258 delete fReco; fReco = 0;
b3ee6797 259 delete fTrClassNamesArr;
296ea9b4 260 delete fSelTracks;
b6c599fe 261 delete fSelPrimTracks;
3a952328 262 delete [] fAmpInTrigger;
d595acbb 263 delete [] fHColuRow;
264 delete [] fHColuRowE;
265 delete [] fHCellMult;
296ea9b4 266 delete [] fHCellFreqNoCut;
2e4d8148 267 delete [] fHCellFreqCut100M;
268 delete [] fHCellFreqCut300M;
ea3fd2d5 269}
717fe7de 270
ea3fd2d5 271//________________________________________________________________________
272void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
273{
717fe7de 274 // Create user objects here.
ea3fd2d5 275
f3582e89 276 cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
277 cout << " fCentVar: " << fCentVar << endl;
278 cout << " fCentFrom: " << fCentFrom << endl;
279 cout << " fCentTo: " << fCentTo << endl;
280 cout << " fVtxZMin: " << fVtxZMin << endl;
281 cout << " fVtxZMax: " << fVtxZMax << endl;
282 cout << " fUseQualFlag: " << fUseQualFlag << endl;
283 cout << " fClusName: \"" << fClusName << "\"" << endl;
284 cout << " fDoNtuple: " << fDoNtuple << endl;
285 cout << " fDoAfterburner: " << fDoAfterburner << endl;
286 cout << " fAsymMax: " << fAsymMax << endl;
287 cout << " fNminCells: " << fNminCells << endl;
288 cout << " fMinE: " << fMinE << endl;
289 cout << " fMinErat: " << fMinErat << endl;
290 cout << " fMinEcc: " << fMinEcc << endl;
291 cout << " fGeoName: \"" << fGeoName << "\"" << endl;
292 cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
38727e64 293 cout << " fIsoDist: " << fIsoDist << endl;
f3582e89 294 cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
f3582e89 295 cout << " fTrCuts: " << fTrCuts << endl;
296 cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
297 cout << " fDoTrMatGeom: " << fDoTrMatGeom << endl;
298 cout << " fTrainMode: " << fTrainMode << endl;
299 cout << " fMarkCells: " << fMarkCells << endl;
300 cout << " fMinL0Time: " << fMinL0Time << endl;
301 cout << " fMaxL0Time: " << fMaxL0Time << endl;
38727e64 302 cout << " fMcMode: " << fMcMode << endl;
cfd7d5b2 303 cout << " fEmbedMode: " << fEmbedMode << endl;
38727e64 304 cout << " fGeom: " << fGeom << endl;
305 cout << " fReco: " << fReco << endl;
807016ea 306 cout << " fDoPSel: " << fDoPSel << endl;
38727e64 307
308 if (!fGeom)
80d69b70 309 fGeom = AliEMCALGeometry::GetInstance(fGeoName);
38727e64 310 else {
311 if (fGeom->GetMatrixForSuperModule(0))
312 fIsGeoMatsSet = kTRUE;
313 }
314 if (!fReco)
315 fReco = new AliEMCALRecoUtils();
b3ee6797 316 fTrClassNamesArr = fTrClassNames.Tokenize(" ");
717fe7de 317 fOutput = new TList();
4ea96211 318 fOutput->SetOwner(1);
296ea9b4 319 fSelTracks = new TObjArray;
b6c599fe 320 fSelPrimTracks = new TObjArray;
ea3fd2d5 321
f5d4ab70 322 if (fDoNtuple) {
a9c1389d 323 TFile *f = OpenFile(1);
324 TDirectory::TContext context(f);
6bf90832 325 if (f) {
326 f->SetCompressionLevel(2);
788ca675 327 fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
6bf90832 328 fNtuple->SetDirectory(f);
3a952328 329 if (fTrainMode) {
330 fNtuple->SetAutoFlush(-2*1024*1024);
331 fNtuple->SetAutoSave(0);
332 } else {
333 fNtuple->SetAutoFlush(-32*1024*1024);
334 fNtuple->SetAutoSave(0);
335 }
336
788ca675 337 fHeader = new AliStaHeader;
338 fNtuple->Branch("header", &fHeader, 16*1024, 99);
339 fPrimVert = new AliStaVertex;
340 fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
341 fSpdVert = new AliStaVertex;
342 fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
343 fTpcVert = new AliStaVertex;
344 fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
345 if (TClass::GetClass("AliStaCluster"))
346 TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
3a952328 347 fClusters = new TClonesArray("AliStaCluster");
788ca675 348 fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
3a952328 349 if (TClass::GetClass("AliStaTrigger"))
350 TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
351 fTriggers = new TClonesArray("AliStaTrigger");
352 fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
cfd7d5b2 353 if (fMcMode||fEmbedMode) {
807016ea 354 if (TClass::GetClass("AliStaPart"))
355 TClass::GetClass("AliStaPart")->IgnoreTObjectStreamer();
356 fMcParts = new TClonesArray("AliStaPart");
357 fNtuple->Branch("mcparts", &fMcParts, 8*16*1024, 99);
358 }
6bf90832 359 }
f5d4ab70 360 }
361
002dcebe 362 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 363 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
364 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
002dcebe 365
d595acbb 366 // histograms
587ce6c6 367 Bool_t th1 = TH1::GetDefaultSumw2();
a49742b5 368 TH1::SetDefaultSumw2(kTRUE);
587ce6c6 369 Bool_t th2 = TH2::GetDefaultSumw2();
a49742b5 370 TH2::SetDefaultSumw2(kTRUE);
b3ee6797 371 fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
372 fHCuts->GetXaxis()->SetBinLabel(1,"All");
373 fHCuts->GetXaxis()->SetBinLabel(2,"PS");
374 fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
375 fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
376 fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
fa443410 377 fOutput->Add(fHCuts);
378 fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
379 fHVertexZ->SetXTitle("z [cm]");
380 fOutput->Add(fHVertexZ);
76332037 381 fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
382 fHVertexZ2->SetXTitle("z [cm]");
383 fOutput->Add(fHVertexZ2);
f2b8fca6 384 fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
fa443410 385 fHCent->SetXTitle(fCentVar.Data());
76332037 386 fOutput->Add(fHCent);
f2b8fca6 387 fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
fa443410 388 fHCentQual->SetXTitle(fCentVar.Data());
389 fOutput->Add(fHCentQual);
2e4d8148 390 fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
391 fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
b3ee6797 392 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
393 const char *name = fTrClassNamesArr->At(i)->GetName();
394 fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
395 fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
396 }
397 fOutput->Add(fHTclsBeforeCuts);
398 fOutput->Add(fHTclsAfterCuts);
90d5b88b 399
d595acbb 400 // histograms for cells
296ea9b4 401 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
402 fHColuRow = new TH2*[nsm];
403 fHColuRowE = new TH2*[nsm];
404 fHCellMult = new TH1*[nsm];
405 for (Int_t i = 0; i<nsm; ++i) {
2e4d8148 406 fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
d595acbb 407 fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
408 fHColuRow[i]->SetXTitle("col (i#eta)");
409 fHColuRow[i]->SetYTitle("row (i#phi)");
2e4d8148 410 fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
90d5b88b 411 fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
d595acbb 412 fHColuRowE[i]->SetXTitle("col (i#eta)");
413 fHColuRowE[i]->SetYTitle("row (i#phi)");
414 fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
415 fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
90d5b88b 416 fHCellMult[i]->SetXTitle("# of cells");
d595acbb 417 fOutput->Add(fHColuRow[i]);
418 fOutput->Add(fHColuRowE[i]);
419 fOutput->Add(fHCellMult[i]);
420 }
2e4d8148 421 fHCellE = new TH1F("hCellE","",250,0.,25.);
d595acbb 422 fHCellE->SetXTitle("E_{cell} [GeV]");
423 fOutput->Add(fHCellE);
fb5a0b69 424 fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
6eb6260e 425 fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
d595acbb 426 fOutput->Add(fHCellH);
fb5a0b69 427 fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
6eb6260e 428 fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
429 fOutput->Add(fHCellM);
fb5a0b69 430 fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
f4ec884e 431 fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
a49742b5 432 fOutput->Add(fHCellM2);
90d5b88b 433
2e4d8148 434 fHCellFreqNoCut = new TH1*[nsm];
435 fHCellFreqCut100M = new TH1*[nsm];
436 fHCellFreqCut300M = new TH1*[nsm];
437 fHCellFreqE = new TH1*[nsm];
296ea9b4 438 for (Int_t i = 0; i<nsm; ++i){
439 Double_t lbin = i*24*48-0.5;
440 Double_t hbin = lbin+24*48;
2e4d8148 441 fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i),
442 Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin);
443 fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
444 Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
445 fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
446 Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
447 fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i),
448 Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
296ea9b4 449 fOutput->Add(fHCellFreqNoCut[i]);
2e4d8148 450 fOutput->Add(fHCellFreqCut100M[i]);
451 fOutput->Add(fHCellFreqCut300M[i]);
452 fOutput->Add(fHCellFreqE[i]);
296ea9b4 453 }
3a952328 454 if (!fMarkCells.IsNull()) {
296ea9b4 455 fHCellCheckE = new TH1*[24*48*nsm];
408dc04c 456 memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
3a952328 457 TObjArray *cells = fMarkCells.Tokenize(" ");
458 Int_t n = cells->GetEntries();
459 Int_t *tcs = new Int_t[n];
460 for (Int_t i=0;i<n;++i) {
461 TString name(cells->At(i)->GetName());
462 tcs[i]=name.Atoi();
463 }
464 for (Int_t i = 0; i<n; ++i) {
296ea9b4 465 Int_t c=tcs[i];
466 if (c<24*48*nsm) {
3a952328 467 fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
296ea9b4 468 fOutput->Add(fHCellCheckE[i]);
469 }
470 }
3a952328 471 delete cells;
472 delete [] tcs;
296ea9b4 473 }
474
d595acbb 475 // histograms for clusters
3a952328 476 if (!fTrainMode) {
477 fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
478 fHClustEccentricity->SetXTitle("#epsilon_{C}");
479 fOutput->Add(fHClustEccentricity);
480 fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
481 fHClustEtaPhi->SetXTitle("#eta");
482 fHClustEtaPhi->SetYTitle("#varphi");
483 fOutput->Add(fHClustEtaPhi);
484 fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
485 fHClustEnergyPt->SetXTitle("E [GeV]");
486 fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
487 fOutput->Add(fHClustEnergyPt);
488 fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
489 fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
490 fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
491 fOutput->Add(fHClustEnergySigma);
492 fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
493 fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
494 fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
495 fOutput->Add(fHClustSigmaSigma);
496 fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
497 fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
498 fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
499 fOutput->Add(fHClustNCellEnergyRatio);
f5e0f1e2 500 fHClustEnergyNCell = new TH2F("hClustEnergyNCell","",200,0,100,50,0,50);
501 fHClustEnergyNCell->SetXTitle("E_{clus}");
502 fHClustEnergyNCell->SetYTitle("N_{cells}");
503 fOutput->Add(fHClustEnergyNCell);
3a952328 504 }
90d5b88b 505
f5e0f1e2 506 // histograms for primary tracks
507 fHPrimTrackPt = new TH1F("hPrimTrackPt",";p_{T} [GeV/c]",500,0,50);
508 fOutput->Add(fHPrimTrackPt);
509 fHPrimTrackEta = new TH1F("hPrimTrackEta",";#eta",40,-2,2);
510 fOutput->Add(fHPrimTrackEta);
511 fHPrimTrackPhi = new TH1F("hPrimTrackPhi",";#varPhi [rad]",63,0,6.3);
512 fOutput->Add(fHPrimTrackPhi);
b6c599fe 513 // histograms for track matching
514 fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
515 fOutput->Add(fHMatchDr);
516 fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
517 fOutput->Add(fHMatchDz);
518 fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
519 fOutput->Add(fHMatchEp);
3a952328 520
d595acbb 521 // histograms for pion candidates
3a952328 522 if (!fTrainMode) {
523 fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
524 fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
525 fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
526 fOutput->Add(fHPionEtaPhi);
527 fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
528 fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
529 fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
530 fOutput->Add(fHPionMggPt);
531 fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
532 fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
533 fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
534 fOutput->Add(fHPionMggAsym);
535 fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
536 fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
537 fHPionMggDgg->SetYTitle("opening angle [grad]");
538 fOutput->Add(fHPionMggDgg);
539 const Int_t nbins = 20;
540 Double_t xbins[nbins] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12.5,15,20,25,50};
541 fPtRanges = new TAxis(nbins-1,xbins);
542 for (Int_t i = 0; i<=nbins; ++i) {
543 fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
544 fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
545 if (i==0)
546 fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
547 else if (i==nbins)
548 fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
549 else
550 fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
551 fOutput->Add(fHPionInvMasses[i]);
552 }
fa443410 553 }
296ea9b4 554
587ce6c6 555 TH1::SetDefaultSumw2(th1);
556 TH2::SetDefaultSumw2(th2);
717fe7de 557 PostData(1, fOutput);
ea3fd2d5 558}
559
560//________________________________________________________________________
561void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
562{
717fe7de 563 // Called for each event.
564
565 if (!InputEvent())
566 return;
567
43cfaa06 568 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
569 fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
44310bce 570 UInt_t offtrigger = 0;
43cfaa06 571 if (fEsdEv) {
572 am->LoadBranch("AliESDRun.");
573 am->LoadBranch("AliESDHeader.");
9809ed8c 574 UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
575 UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
576 Bool_t desc1 = (mask1 >> 18) & 0x1;
577 Bool_t desc2 = (mask2 >> 18) & 0x1;
2cee9afd 578 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(18)=="EMCAL"
9809ed8c 579 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
580 mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
581 mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
44310bce 582 return;
583 }
584 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
43cfaa06 585 } else {
586 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
408dc04c 587 if (!fAodEv) {
588 AliFatal("Neither ESD nor AOD event found");
589 return;
590 }
43cfaa06 591 am->LoadBranch("header");
44310bce 592 offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
593 }
38727e64 594 if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
44310bce 595 AliWarning(Form("EMCAL not in fast only partition"));
596 return;
43cfaa06 597 }
b6c599fe 598 if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
27c2e3d9 599 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
600 AliCDBManager *cdb = AliCDBManager::Instance();
601 if (!cdb->IsDefaultStorageSet())
602 cdb->SetDefaultStorage("raw://");
603 Int_t runno = InputEvent()->GetRunNumber();
604 if (runno != cdb->GetRun())
605 cdb->SetRun(runno);
606 AliGeomManager::LoadGeometry();
607 }
608
609 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
610 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
9809ed8c 611 for (Int_t i=0; i<nsm; ++i) {
612 const TGeoHMatrix *geom = 0;
613 if (fEsdEv)
614 geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
615 else
616 geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
617 if (!geom)
618 continue;
619 geom->Print();
620 fGeom->SetMisalMatrix(geom,i);
27c2e3d9 621 }
622 fIsGeoMatsSet = kTRUE;
623 }
624
625 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
626 if (fEsdEv) {
627 const AliESDRun *erun = fEsdEv->GetESDRun();
628 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
629 erun->GetCurrentDip(),
630 AliMagF::kConvLHC,
631 kFALSE,
632 erun->GetBeamEnergy(),
633 erun->GetBeamType());
634 TGeoGlobalMagField::Instance()->SetField(field);
635 } else {
636 Double_t pol = -1; //polarity
637 Double_t be = -1; //beam energy
638 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
639 Int_t runno = fAodEv->GetRunNumber();
640 if (runno>=136851 && runno<138275) {
641 pol = -1;
642 be = 2760;
643 btype = AliMagF::kBeamTypeAA;
644 } else if (runno>=138275 && runno<=139517) {
645 pol = +1;
646 be = 2760;
647 btype = AliMagF::kBeamTypeAA;
648 } else {
649 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
650 }
651 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
652 }
653 }
654
b3ee6797 655 Int_t cut = 1;
656 fHCuts->Fill(cut++);
657
658 TString trgclasses;
659 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
660 if (h) {
661 trgclasses = fEsdEv->GetFiredTriggerClasses();
662 } else {
663 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
664 if (h2)
665 trgclasses = h2->GetFiredTriggerClasses();
666 }
667 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
668 const char *name = fTrClassNamesArr->At(i)->GetName();
669 if (trgclasses.Contains(name))
670 fHTclsBeforeCuts->Fill(1+i);
671 }
672
807016ea 673 if (fDoPSel && offtrigger==0)
b3ee6797 674 return;
675
fa443410 676 fHCuts->Fill(cut++);
286b47a5 677
717fe7de 678 const AliCentrality *centP = InputEvent()->GetCentrality();
679 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
fa443410 680 fHCent->Fill(cent);
717fe7de 681 if (cent<fCentFrom||cent>fCentTo)
286b47a5 682 return;
43cfaa06 683
fa443410 684 fHCuts->Fill(cut++);
286b47a5 685
43cfaa06 686 if (fUseQualFlag) {
687 if (centP->GetQuality()>0)
688 return;
689 }
286b47a5 690
fa443410 691 fHCentQual->Fill(cent);
692 fHCuts->Fill(cut++);
717fe7de 693
43cfaa06 694 if (fEsdEv) {
fa443410 695 am->LoadBranch("PrimaryVertex.");
696 am->LoadBranch("SPDVertex.");
697 am->LoadBranch("TPCVertex.");
43cfaa06 698 } else {
699 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
700 am->LoadBranch("vertices");
3c661da5 701 if (!fAodEv) return;
43cfaa06 702 }
703
704 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
717fe7de 705 if (!vertex)
706 return;
707
fa443410 708 fHVertexZ->Fill(vertex->GetZ());
286b47a5 709
717fe7de 710 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
711 return;
286b47a5 712
fa443410 713 fHCuts->Fill(cut++);
76332037 714 fHVertexZ2->Fill(vertex->GetZ());
717fe7de 715
b3ee6797 716 // count number of accepted events
717 ++fNEvs;
718
719 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
720 const char *name = fTrClassNamesArr->At(i)->GetName();
721 if (trgclasses.Contains(name))
722 fHTclsAfterCuts->Fill(1+i);
723 }
724
43cfaa06 725 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
c2fe5f0e 726 fDigits = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
b6c599fe 727 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 728 fEsdCells = 0; // will be set if ESD input used
b6c599fe 729 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 730 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
731 fAodCells = 0; // will be set if AOD input used
717fe7de 732
43cfaa06 733 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
734 Bool_t clusattached = 0;
735 Bool_t recalibrated = 0;
736 if (1 && !fClusName.IsNull()) {
737 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
738 TObjArray *ts = am->GetTasks();
739 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
740 if (cltask && cltask->GetClusters()) {
c2fe5f0e 741 fRecPoints = cltask->GetClusters();
742 fDigits = cltask->GetDigits();
43cfaa06 743 clusattached = cltask->GetAttachClusters();
744 if (cltask->GetCalibData()!=0)
745 recalibrated = kTRUE;
746 }
747 }
38727e64 748 if (1 && !fClusName.IsNull()) {
749 TList *l = 0;
750 if (AODEvent())
751 l = AODEvent()->GetList();
752 else if (fAodEv)
753 l = fAodEv->GetList();
717fe7de 754 if (l) {
38727e64 755 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
717fe7de 756 }
ea3fd2d5 757 }
717fe7de 758
759 if (fEsdEv) { // ESD input mode
43cfaa06 760 if (1 && (!fRecPoints||clusattached)) {
761 if (!clusattached)
762 am->LoadBranch("CaloClusters");
717fe7de 763 TList *l = fEsdEv->GetList();
717fe7de 764 if (l) {
38727e64 765 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
ea3fd2d5 766 }
767 }
43cfaa06 768 if (1) {
769 if (!recalibrated)
770 am->LoadBranch("EMCALCells.");
717fe7de 771 fEsdCells = fEsdEv->GetEMCALCells();
772 }
773 } else if (fAodEv) { // AOD input mode
43cfaa06 774 if (1 && (!fAodClusters || clusattached)) {
775 if (!clusattached)
776 am->LoadBranch("caloClusters");
717fe7de 777 TList *l = fAodEv->GetList();
717fe7de 778 if (l) {
38727e64 779 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
ea3fd2d5 780 }
781 }
43cfaa06 782 if (1) {
783 if (!recalibrated)
784 am->LoadBranch("emcalCells");
717fe7de 785 fAodCells = fAodEv->GetEMCALCells();
786 }
787 } else {
788 AliFatal("Impossible to not have either pointer to ESD or AOD event");
789 }
ea3fd2d5 790
43cfaa06 791 if (1) {
792 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
c2fe5f0e 793 AliDebug(2,Form("fDigits set: %p", fDigits));
43cfaa06 794 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
795 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
796 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
797 AliDebug(2,Form("fAodCells set: %p", fAodCells));
798 }
799
a49742b5 800 if (fDoAfterburner)
801 ClusterAfterburner();
6eb6260e 802
38727e64 803 CalcMcInfo();
3a952328 804 CalcCaloTriggers();
b6c599fe 805 CalcPrimTracks();
3a952328 806 CalcTracks();
296ea9b4 807 CalcClusterProps();
808
286b47a5 809 FillCellHists();
3a952328 810 if (!fTrainMode) {
811 FillClusHists();
812 FillPionHists();
813 FillOtherHists();
814 }
38727e64 815 FillMcHists();
788ca675 816 FillNtuple();
ea3fd2d5 817
3a952328 818 if (fTrainMode) {
819 fSelTracks->Clear();
820 fSelPrimTracks->Clear();
cfd7d5b2 821 if (fMcParts)
807016ea 822 fMcParts->Clear();
4ea96211 823 if (fTriggers)
824 fTriggers->Clear();
825 if (fClusters)
826 fClusters->Clear();
3a952328 827 }
de4cee41 828
717fe7de 829 PostData(1, fOutput);
ea3fd2d5 830}
831
832//________________________________________________________________________
833void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
834{
717fe7de 835 // Terminate called at the end of analysis.
f5d4ab70 836
837 if (fNtuple) {
4ea96211 838 TFile *f = OpenFile(1);
a9c1389d 839 TDirectory::TContext context(f);
f5d4ab70 840 if (f)
841 fNtuple->Write();
842 }
d9f26424 843
bc8a45bd 844 AliInfo(Form("%s: Accepted %lld events ", GetName(), fNEvs));
ea3fd2d5 845}
ea3fd2d5 846
296ea9b4 847//________________________________________________________________________
3a952328 848void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
296ea9b4 849{
3a952328 850 // Calculate triggers
296ea9b4 851
38727e64 852 if (fAodEv)
853 return; // information not available in AOD
854
4ea96211 855 if (!fTriggers)
856 return;
857
3a952328 858 fTriggers->Clear();
296ea9b4 859
3a952328 860 AliVCaloCells *cells = fEsdCells;
861 if (!cells)
862 cells = fAodCells;
863 if (!cells)
864 return;
865
866 Int_t ncells = cells->GetNumberOfCells();
867 if (ncells<=0)
868 return;
869
870 if (ncells>fNAmpInTrigger) {
871 delete [] fAmpInTrigger;
872 fAmpInTrigger = new Float_t[ncells];
873 fNAmpInTrigger = ncells;
296ea9b4 874 }
3a952328 875 for (Int_t i=0;i<ncells;++i)
876 fAmpInTrigger[i] = 0;
296ea9b4 877
3a952328 878 std::map<Short_t,Short_t> map;
879 for (Short_t pos=0;pos<ncells;++pos) {
880 Short_t id = cells->GetCellNumber(pos);
881 map[id]=pos;
882 }
883
884 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
885 am->LoadBranch("EMCALTrigger.");
886
887 AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
888 if (!triggers)
889 return;
890 if (triggers->GetEntries()<=0)
b6c599fe 891 return;
892
3a952328 893 triggers->Reset();
894 Int_t ntrigs=0;
895 while (triggers->Next()) {
896 Int_t gCol=0, gRow=0, ntimes=0;
897 triggers->GetPosition(gCol,gRow);
898 triggers->GetNL0Times(ntimes);
899 if (ntimes<1)
900 continue;
901 Float_t amp=0;
902 triggers->GetAmplitude(amp);
903 Int_t find = -1;
904 fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
905 if (find<0)
906 continue;
907 Int_t cidx[4] = {-1};
908 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
909 if (!ret)
910 continue;
911 Int_t trgtimes[25];
912 triggers->GetL0Times(trgtimes);
913 Int_t mintime = trgtimes[0];
914 Int_t maxtime = trgtimes[0];
915 Bool_t trigInTimeWindow = 0;
916 for (Int_t i=0;i<ntimes;++i) {
917 if (trgtimes[i]<mintime)
918 mintime = trgtimes[i];
919 if (maxtime<trgtimes[i])
920 maxtime = trgtimes[i];
f2b8fca6 921 if ((fMinL0Time<=trgtimes[i]) && (fMaxL0Time>=trgtimes[i]))
3a952328 922 trigInTimeWindow = 1;
923 }
924
925 Double_t tenergy = 0;
926 Double_t tphi=0;
927 Double_t teta=0;
928 for (Int_t i=0;i<3;++i) {
929 Short_t pos = -1;
930 std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
931 if (it!=map.end())
932 pos = it->second;
933 if (pos<0)
b6c599fe 934 continue;
3a952328 935 if (trigInTimeWindow)
5fe1ca23 936 fAmpInTrigger[pos] = amp;
3a952328 937 Float_t eta=-1, phi=-1;
938 fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
939 Double_t en= cells->GetAmplitude(pos);
940 tenergy+=en;
941 teta+=eta*en;
942 tphi+=phi*en;
943 }
c47674e0 944
945 if (tenergy<=0)
946 continue;
947
3a952328 948 teta/=tenergy;
949 tphi/=tenergy;
c47674e0 950
3a952328 951 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
952 trignew->fE = tenergy;
953 trignew->fEta = teta;
954 trignew->fPhi = tphi;
955 trignew->fAmp = amp;
956 trignew->fMinTime = mintime;
957 trignew->fMaxTime = maxtime;
958 }
959}
960
961//________________________________________________________________________
962void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
963{
964 // Calculate cluster properties
965
4ea96211 966 if (!fClusters)
967 return;
968
3a952328 969 fClusters->Clear();
970
971 AliVCaloCells *cells = fEsdCells;
972 if (!cells)
973 cells = fAodCells;
974 if (!cells)
975 return;
976
977 TObjArray *clusters = fEsdClusters;
978 if (!clusters)
979 clusters = fAodClusters;
980 if (!clusters)
981 return;
982
983 Int_t ncells = cells->GetNumberOfCells();
984 Int_t nclus = clusters->GetEntries();
985 Int_t ntrks = fSelTracks->GetEntries();
587ce6c6 986 Int_t btracks[6][ntrks];
987 memset(btracks,0,sizeof(Int_t)*6*ntrks);
3a952328 988
989 std::map<Short_t,Short_t> map;
990 for (Short_t pos=0;pos<ncells;++pos) {
991 Short_t id = cells->GetCellNumber(pos);
992 map[id]=pos;
993 }
994
8c56d760 995 TObjArray filtMcParts;
cfd7d5b2 996 if (fMcParts) {
8c56d760 997 Int_t nmc = fMcParts->GetEntries();
998 for (Int_t i=0; i<nmc; ++i) {
999 AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
1000 if (pa->OnEmcal())
1001 filtMcParts.Add(pa);
1002 }
1003 }
1004
3a952328 1005 for(Int_t i=0, ncl=0; i<nclus; ++i) {
1006 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
b5ba4932 1007
3a952328 1008 if (!clus)
1009 continue;
1010 if (!clus->IsEMCAL())
1011 continue;
1012 if (clus->E()<fMinE)
1013 continue;
1014
1015 Float_t clsPos[3] = {0};
1016 clus->GetPosition(clsPos);
1017 TVector3 clsVec(clsPos);
1018 Double_t vertex[3] = {0};
1019 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1020 TLorentzVector clusterVec;
1021 clus->GetMomentum(clusterVec,vertex);
1022 Double_t clsEta = clusterVec.Eta();
1023
1024 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
1025 cl->fE = clus->E();
1026 cl->fR = clsVec.Perp();
1027 cl->fEta = clsVec.Eta();
1028 cl->fPhi = clsVec.Phi();
1029 cl->fN = clus->GetNCells();
1030 cl->fN1 = GetNCells(clus,0.100);
1031 cl->fN3 = GetNCells(clus,0.300);
1032 Short_t id = -1;
1033 Double_t emax = GetMaxCellEnergy(clus, id);
1034 cl->fIdMax = id;
1035 cl->fEmax = emax;
c2fe5f0e 1036 cl->fTmax = cells->GetCellTime(id);
3a952328 1037 if (clus->GetDistanceToBadChannel()<10000)
1038 cl->fDbc = clus->GetDistanceToBadChannel();
1039 if (!TMath::IsNaN(clus->GetDispersion()))
1040 cl->fDisp = clus->GetDispersion();
1041 if (!TMath::IsNaN(clus->GetM20()))
0fbe8d4f 1042 cl->fM20 = clus->GetM20();
3a952328 1043 if (!TMath::IsNaN(clus->GetM02()))
0fbe8d4f 1044 cl->fM02 = clus->GetM02();
3a952328 1045 Double_t maxAxis = 0, minAxis = 0;
1046 GetSigma(clus,maxAxis,minAxis);
1047 clus->SetTOF(maxAxis); // store sigma in TOF
1048 cl->fSig = maxAxis;
1049 Double_t clusterEcc = 0;
1050 if (maxAxis > 0)
1051 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
1052 clus->SetChi2(clusterEcc); // store ecc in chi2
1053 cl->fEcc = clusterEcc;
1054 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
1055 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
1056 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
1057 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
1058 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
0fbe8d4f 1059
38727e64 1060 if (fAmpInTrigger) { // fill trigger info if present
1061 Double_t trigpen = 0;
1062 Double_t trignen = 0;
1063 for(Int_t j=0; j<cl->fN; ++j) {
1064 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1065 Short_t pos = -1;
1066 std::map<Short_t,Short_t>::iterator it = map.find(cid);
1067 if (it!=map.end())
1068 pos = it->second;
1069 if (pos<0)
1070 continue;
1071 if (fAmpInTrigger[pos]>0)
1072 trigpen += cells->GetAmplitude(pos);
1073 else if (fAmpInTrigger[pos]<0)
1074 trignen += cells->GetAmplitude(pos);
1075 }
1076 if (trigpen>0) {
1077 cl->fIsTrigM = 1;
1078 cl->fTrigE = trigpen;
1079 }
1080 if (trignen>0) {
1081 cl->fIsTrigM = 1;
1082 cl->fTrigMaskE = trignen;
1083 }
3a952328 1084 }
0fbe8d4f 1085 cl->fIsShared = IsShared(clus);
3a952328 1086
1087 // track matching
1088 Double_t mind2 = 1e10;
1089 for(Int_t j = 0; j<ntrks; ++j) {
1090 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
b6c599fe 1091 if (!track)
1092 continue;
3a952328 1093
1094 if (TMath::Abs(clsEta-track->Eta())>0.5)
b6c599fe 1095 continue;
3a952328 1096
1097 TVector3 vec(clsPos);
1098 Int_t index = (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
1099 if (btracks[index-4][j]) {
b6c599fe 1100 continue;
3a952328 1101 }
b6c599fe 1102
3a952328 1103 Float_t tmpR=-1, tmpZ=-1;
1104 if (!fDoTrMatGeom) {
1105 AliExternalTrackParam *tParam = 0;
1106 if (fEsdEv) {
1107 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
1108 tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
1109 } else
1110 tParam = new AliExternalTrackParam(track);
1111
1112 Double_t bfield[3] = {0};
1113 track->GetBxByBz(bfield);
1114 Double_t alpha = (index+0.5)*20*TMath::DegToRad();
1115 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1116 tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1117 Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
1118 if (!ret) {
1119 btracks[index-4][j]=1;
1120 delete tParam;
1121 continue;
b6c599fe 1122 }
3a952328 1123 Double_t trkPos[3] = {0};
1124 tParam->GetXYZ(trkPos); //Get the extrapolated global position
38727e64 1125 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2) +
1126 TMath::Power(clsPos[1]-trkPos[1],2) +
1127 TMath::Power(clsPos[2]-trkPos[2],2) );
3a952328 1128 tmpZ = clsPos[2]-trkPos[2];
1129 delete tParam;
1130 } else {
1131 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
1132 continue;
1133 AliExternalTrackParam tParam(track);
1134 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
1135 continue;
b6c599fe 1136 }
3a952328 1137
1138 Double_t d2 = tmpR;
1139 if (mind2>d2) {
1140 mind2=d2;
1141 cl->fTrDz = tmpZ;
1142 cl->fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
1143 cl->fTrEp = clus->E()/track->P();
f3582e89 1144 cl->fIsTrackM = 1;
3a952328 1145 }
1146 }
1147
f3582e89 1148 if (cl->fIsTrackM) {
3a952328 1149 fHMatchDr->Fill(cl->fTrDr);
1150 fHMatchDz->Fill(cl->fTrDz);
1151 fHMatchEp->Fill(cl->fTrEp);
b6c599fe 1152 }
8c56d760 1153
1154 //mc matching
cfd7d5b2 1155 if (fMcParts) {
8c56d760 1156 Int_t nmc = filtMcParts.GetEntries();
1157 Double_t diffR2 = 1e9;
1158 AliStaPart *msta = 0;
1159 for (Int_t j=0; j<nmc; ++j) {
1160 AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
1161 Double_t t1=clsVec.Eta()-pa->fVEta;
1162 Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
1163 Double_t tmp = t1*t1+t2*t2;
1164 if (tmp<diffR2) {
1165 diffR2 = tmp;
1166 msta = pa;
1167 }
1168 }
c2fe5f0e 1169 if (diffR2<10 && msta!=0) {
1170 cl->fMcLabel = msta->fLab;
1171 }
1172 }
1173
1174 cl->fEmbE = 0;
1175 if (fDigits && fEmbedMode) {
1176 for(Int_t j=0; j<cl->fN; ++j) {
1177 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1178 Short_t pos = -1;
1179 std::map<Short_t,Short_t>::iterator it = map.find(cid);
1180 if (it!=map.end())
1181 pos = it->second;
1182 if (pos<0)
1183 continue;
1184 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
1185 if (!digit)
1186 continue;
1187 if (digit->GetId() != cid) {
1188 AliError(Form("Ids should be equal: %d %d", cid, digit->GetId()));
1189 continue;
1190 }
1191 if (digit->GetType()<-1) {
1192 cl->fEmbE += digit->GetChi2();
1193 }
1194 }
8c56d760 1195 }
b6c599fe 1196 }
1197}
1198
1199//________________________________________________________________________
1200void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
1201{
a2de5ca1 1202 // Calculate track properties for primary tracks.
b6c599fe 1203
1204 fSelPrimTracks->Clear();
1205
1206 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1207 TClonesArray *tracks = 0;
1208 if (fEsdEv) {
1209 am->LoadBranch("Tracks");
1210 TList *l = fEsdEv->GetList();
1211 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1212 } else {
1213 am->LoadBranch("tracks");
1214 TList *l = fAodEv->GetList();
1215 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1216 }
1217
296ea9b4 1218 if (!tracks)
1219 return;
1220
0ec74551 1221 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 1222 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1223 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
0ec74551 1224
296ea9b4 1225 if (fEsdEv) {
b6c599fe 1226 fSelPrimTracks->SetOwner(kTRUE);
0ec74551 1227 am->LoadBranch("PrimaryVertex.");
1228 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1229 am->LoadBranch("SPDVertex.");
1230 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1231 am->LoadBranch("Tracks");
1232 const Int_t Ntracks = tracks->GetEntries();
1233 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1234 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1235 if (!track) {
1236 AliWarning(Form("Could not receive track %d\n", iTracks));
1237 continue;
1238 }
1239 if (fTrCuts && !fTrCuts->IsSelected(track))
296ea9b4 1240 continue;
1241 Double_t eta = track->Eta();
1242 if (eta<-1||eta>1)
1243 continue;
0ec74551 1244 if (track->Phi()<phimin||track->Phi()>phimax)
1245 continue;
2e4d8148 1246
0ec74551 1247 AliESDtrack copyt(*track);
1248 Double_t bfield[3];
1249 copyt.GetBxByBz(bfield);
1250 AliExternalTrackParam tParam;
1251 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1252 if (!relate)
1253 continue;
1254 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
2e4d8148 1255
0ec74551 1256 Double_t p[3] = { 0. };
1257 copyt.GetPxPyPz(p);
1258 Double_t pos[3] = { 0. };
1259 copyt.GetXYZ(pos);
1260 Double_t covTr[21] = { 0. };
1261 copyt.GetCovarianceXYZPxPyPz(covTr);
1262 Double_t pid[10] = { 0. };
1263 copyt.GetESDpid(pid);
1264 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1265 copyt.GetLabel(),
1266 p,
1267 kTRUE,
1268 pos,
1269 kFALSE,
1270 covTr,
1271 (Short_t)copyt.GetSign(),
1272 copyt.GetITSClusterMap(),
1273 pid,
1274 0,/*fPrimaryVertex,*/
1275 kTRUE, // check if this is right
1276 vtx->UsesTrack(copyt.GetID()));
1277 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1278 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1279 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1280 if(ndf>0)
1281 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1282 else
1283 aTrack->SetChi2perNDF(-1);
1284 aTrack->SetFlags(copyt.GetStatus());
1285 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
b6c599fe 1286 fSelPrimTracks->Add(aTrack);
296ea9b4 1287 }
1288 } else {
1289 Int_t ntracks = tracks->GetEntries();
1290 for (Int_t i=0; i<ntracks; ++i) {
1291 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1292 if (!track)
1293 continue;
1294 Double_t eta = track->Eta();
1295 if (eta<-1||eta>1)
1296 continue;
0ec74551 1297 if (track->Phi()<phimin||track->Phi()>phimax)
1298 continue;
b6c599fe 1299 if(track->GetTPCNcls()<fMinNClusPerTr)
296ea9b4 1300 continue;
b6c599fe 1301 //todo: Learn how to set/filter AODs for prim/sec tracks
1302 fSelPrimTracks->Add(track);
296ea9b4 1303 }
1304 }
1305}
1306
1307//________________________________________________________________________
3a952328 1308void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
296ea9b4 1309{
3a952328 1310 // Calculate track properties (including secondaries).
b6c599fe 1311
3a952328 1312 fSelTracks->Clear();
296ea9b4 1313
3a952328 1314 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1315 TClonesArray *tracks = 0;
1316 if (fEsdEv) {
1317 am->LoadBranch("Tracks");
1318 TList *l = fEsdEv->GetList();
1319 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1320 } else {
1321 am->LoadBranch("tracks");
1322 TList *l = fAodEv->GetList();
1323 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1324 }
296ea9b4 1325
3a952328 1326 if (!tracks)
1327 return;
296ea9b4 1328
3a952328 1329 if (fEsdEv) {
1330 const Int_t Ntracks = tracks->GetEntries();
1331 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1332 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1333 if (!track) {
1334 AliWarning(Form("Could not receive track %d\n", iTracks));
1335 continue;
1336 }
1337 if (fTrCuts && !fTrCuts->IsSelected(track))
1338 continue;
1339 Double_t eta = track->Eta();
1340 if (eta<-1||eta>1)
1341 continue;
1342 fSelTracks->Add(track);
1343 }
1344 } else {
1345 Int_t ntracks = tracks->GetEntries();
1346 for (Int_t i=0; i<ntracks; ++i) {
1347 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
296ea9b4 1348 if (!track)
1349 continue;
3a952328 1350 Double_t eta = track->Eta();
1351 if (eta<-1||eta>1)
c4236a58 1352 continue;
3a952328 1353 if(track->GetTPCNcls()<fMinNClusPerTr)
b6c599fe 1354 continue;
2e4d8148 1355
3a952328 1356 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
2e4d8148 1357 AliExternalTrackParam tParam(track);
3a952328 1358 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1359 Double_t trkPos[3];
1360 tParam.GetXYZ(trkPos);
1361 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1362 }
b6c599fe 1363 }
3a952328 1364 fSelTracks->Add(track);
c4236a58 1365 }
296ea9b4 1366 }
1367}
1368
323834f0 1369//________________________________________________________________________
3a952328 1370void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
323834f0 1371{
296ea9b4 1372 // Run custer reconstruction afterburner.
323834f0 1373
1374 AliVCaloCells *cells = fEsdCells;
1375 if (!cells)
1376 cells = fAodCells;
1377
1378 if (!cells)
1379 return;
1380
1381 Int_t ncells = cells->GetNumberOfCells();
1382 if (ncells<=0)
1383 return;
1384
1385 Double_t cellMeanE = 0, cellSigE = 0;
1386 for (Int_t i = 0; i<ncells; ++i) {
1387 Double_t cellE = cells->GetAmplitude(i);
1388 cellMeanE += cellE;
1389 cellSigE += cellE*cellE;
1390 }
1391 cellMeanE /= ncells;
1392 cellSigE /= ncells;
1393 cellSigE -= cellMeanE*cellMeanE;
1394 if (cellSigE<0)
1395 cellSigE = 0;
1396 cellSigE = TMath::Sqrt(cellSigE / ncells);
1397
1398 Double_t subE = cellMeanE - 7*cellSigE;
1399 if (subE<0)
1400 return;
1401
1402 for (Int_t i = 0; i<ncells; ++i) {
1403 Short_t id=-1;
1404 Double_t amp=0,time=0;
1405 if (!cells->GetCell(i, id, amp, time))
1406 continue;
1407 amp -= cellMeanE;
1408 if (amp<0.001)
1409 amp = 0;
1410 cells->SetCell(i, id, amp, time);
1411 }
1412
1413 TObjArray *clusters = fEsdClusters;
1414 if (!clusters)
1415 clusters = fAodClusters;
323834f0 1416 if (!clusters)
1417 return;
1418
1419 Int_t nclus = clusters->GetEntries();
1420 for (Int_t i = 0; i<nclus; ++i) {
1421 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1422 if (!clus->IsEMCAL())
1423 continue;
1424 Int_t nc = clus->GetNCells();
1425 Double_t clusE = 0;
1426 UShort_t ids[100] = {0};
1427 Double_t fra[100] = {0};
1428 for (Int_t j = 0; j<nc; ++j) {
1429 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1430 Double_t cen = cells->GetCellAmplitude(id);
1431 clusE += cen;
1432 if (cen>0) {
1433 ids[nc] = id;
1434 ++nc;
1435 }
1436 }
1437 if (clusE<=0) {
1438 clusters->RemoveAt(i);
1439 continue;
1440 }
1441
1442 for (Int_t j = 0; j<nc; ++j) {
1443 Short_t id = ids[j];
1444 Double_t cen = cells->GetCellAmplitude(id);
1445 fra[j] = cen/clusE;
1446 }
1447 clus->SetE(clusE);
1448 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1449 if (aodclus) {
1450 aodclus->Clear("");
1451 aodclus->SetNCells(nc);
1452 aodclus->SetCellsAmplitudeFraction(fra);
1453 aodclus->SetCellsAbsId(ids);
1454 }
1455 }
1456 clusters->Compress();
1457}
1458
286b47a5 1459//________________________________________________________________________
1460void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1461{
1462 // Fill histograms related to cell properties.
d595acbb 1463
90d5b88b 1464 AliVCaloCells *cells = fEsdCells;
1465 if (!cells)
1466 cells = fAodCells;
1467
296ea9b4 1468 if (!cells)
1469 return;
90d5b88b 1470
296ea9b4 1471 Int_t cellModCount[12] = {0};
1472 Double_t cellMaxE = 0;
1473 Double_t cellMeanE = 0;
1474 Int_t ncells = cells->GetNumberOfCells();
1475 if (ncells==0)
1476 return;
1477
1478 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1479
1480 for (Int_t i = 0; i<ncells; ++i) {
1481 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1482 Double_t cellE = cells->GetAmplitude(i);
1483 fHCellE->Fill(cellE);
1484 if (cellE>cellMaxE)
1485 cellMaxE = cellE;
1486 cellMeanE += cellE;
1487
1488 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1489 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1490 if (!ret) {
1491 AliError(Form("Could not get cell index for %d", absID));
1492 continue;
1493 }
1494 ++cellModCount[iSM];
1495 Int_t iPhi=-1, iEta=-1;
1496 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1497 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1498 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1499 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1500 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1501 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1502 if (fHCellCheckE && fHCellCheckE[absID])
1503 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1504 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1505 }
1506 fHCellH->Fill(cellMaxE);
1507 cellMeanE /= ncells;
1508 fHCellM->Fill(cellMeanE);
1509 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1510 for (Int_t i=0; i<nsm; ++i)
1511 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1512}
1513
1514//________________________________________________________________________
1515void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1516{
90d5b88b 1517 // Fill histograms related to cluster properties.
fa443410 1518
90d5b88b 1519 TObjArray *clusters = fEsdClusters;
1520 if (!clusters)
1521 clusters = fAodClusters;
296ea9b4 1522 if (!clusters)
1523 return;
90d5b88b 1524
296ea9b4 1525 Double_t vertex[3] = {0};
1526 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1527
1528 Int_t nclus = clusters->GetEntries();
1529 for(Int_t i = 0; i<nclus; ++i) {
1530 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1531 if (!clus)
1532 continue;
1533 if (!clus->IsEMCAL())
1534 continue;
90d5b88b 1535 TLorentzVector clusterVec;
296ea9b4 1536 clus->GetMomentum(clusterVec,vertex);
3a952328 1537 Double_t maxAxis = clus->GetTOF(); //sigma
1538 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1539 fHClustEccentricity->Fill(clusterEcc);
1540 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1541 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1542 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1543 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1544 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
f5e0f1e2 1545 //====
1546 fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
788ca675 1547 }
1548}
b3ee6797 1549
38727e64 1550//________________________________________________________________________
1551void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1552{
1553 // Get Mc truth particle information.
1554
cfd7d5b2 1555 if (!fMcParts)
38727e64 1556 return;
1557
807016ea 1558 fMcParts->Clear();
1559
1560 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1561 Double_t etamin = -0.7;
1562 Double_t etamax = +0.7;
1563 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1564 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1565
56fd6cb2 1566 if (fAodEv) {
807016ea 1567 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1568 am->LoadBranch(AliAODMCParticle::StdBranchName());
56fd6cb2 1569 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1570 if (!tca)
1571 return;
1572
1573 Int_t nents = tca->GetEntries();
807016ea 1574 for(int it=0; it<nents; ++it) {
56fd6cb2 1575 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
807016ea 1576 part->Print();
56fd6cb2 1577
b5ba4932 1578 // pion or eta meson or direct photon
56fd6cb2 1579 if(part->GetPdgCode() == 111) {
1580 } else if(part->GetPdgCode() == 221) {
b5ba4932 1581 } else if(part->GetPdgCode() == 22 ) {
1582 } else
56fd6cb2 1583 continue;
1584
1585 // primary particle
1586 Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1587 if(dR > 1.0)
1588 continue;
56fd6cb2 1589
807016ea 1590 // kinematic cuts
1591 Double_t pt = part->Pt() ;
1592 if (pt<0.5)
1593 continue;
1594 Double_t eta = part->Eta();
1595 if (eta<etamin||eta>etamax)
1596 continue;
1597 Double_t phi = part->Phi();
1598 if (phi<phimin||phi>phimax)
1599 continue;
1600
1601 ProcessDaughters(part, it, tca);
56fd6cb2 1602 }
807016ea 1603 return;
1604 }
38727e64 1605
1606 AliMCEvent *mcEvent = MCEvent();
1607 if (!mcEvent)
1608 return;
1609
1610 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1611 if (!evtVtx)
1612 return;
1613
1614 mcEvent->PreReadAll();
38727e64 1615
1616 Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1617 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1618 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1619 if (!mcP)
1620 continue;
1621
b5ba4932 1622 // pion or eta meson or direct photon
807016ea 1623 if(mcP->PdgCode() == 111) {
1624 } else if(mcP->PdgCode() == 221) {
b5ba4932 1625 } else if(mcP->PdgCode() == 22 ) {
38727e64 1626 } else
1627 continue;
1628
1629 // primary particle
807016ea 1630 Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1631 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
38727e64 1632 if(dR > 1.0)
1633 continue;
1634
38727e64 1635 // kinematic cuts
807016ea 1636 Double_t pt = mcP->Pt() ;
1637 if (pt<0.5)
38727e64 1638 continue;
807016ea 1639 Double_t eta = mcP->Eta();
1640 if (eta<etamin||eta>etamax)
38727e64 1641 continue;
807016ea 1642 Double_t phi = mcP->Phi();
1643 if (phi<phimin||phi>phimax)
1644 continue;
1645
1646 ProcessDaughters(mcP, iTrack, mcEvent);
38727e64 1647 }
1648}
1649
788ca675 1650//________________________________________________________________________
1651void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1652{
1653 // Fill ntuple.
1654
1655 if (!fNtuple)
1656 return;
1657
1658 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1659 if (fAodEv) {
1660 fHeader->fRun = fAodEv->GetRunNumber();
1661 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1662 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1663 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1664 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1665 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1666 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1667 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1668 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1669 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1670 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1671 } else {
1672 fHeader->fRun = fEsdEv->GetRunNumber();
1673 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1674 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1675 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1676 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1677 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1678 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1679 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1680 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1681 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1682 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1683 Float_t v0CorrR = 0;
1684 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1685 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1686 if (mult)
1687 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1688 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
c0d2e1f2 1689 AliTriggerAnalysis trAn; /// Trigger Analysis
1690 Bool_t v0B = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0C);
1691 Bool_t v0A = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0A);
1692 fHeader->fV0And = v0A && v0B;
788ca675 1693 }
1694 AliCentrality *cent = InputEvent()->GetCentrality();
1695 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1696 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1697 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1698 fHeader->fCqual = cent->GetQuality();
b6c599fe 1699
1700 AliEventplane *ep = InputEvent()->GetEventplane();
1701 if (ep) {
1702 if (ep->GetQVector())
1703 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1704 else
1705 fHeader->fPsi = -1;
b6c599fe 1706 if (ep->GetQsub1()&&ep->GetQsub2())
1707 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1708 else
1709 fHeader->fPsiRes = 0;
1710 }
1711
788ca675 1712 Double_t val = 0;
1713 TString trgclasses(fHeader->fFiredTriggers);
1714 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1715 const char *name = fTrClassNamesArr->At(j)->GetName();
1716 if (trgclasses.Contains(name))
1717 val += TMath::Power(2,j);
1718 }
1719 fHeader->fTcls = (UInt_t)val;
1720
5fe1ca23 1721 fHeader->fNSelTr = fSelTracks->GetEntries();
1722 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
f5e0f1e2 1723 fHeader->fNSelPrimTr1 = 0;
1724 fHeader->fNSelPrimTr2 = 0;
1725 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1726 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1727 if(track->Pt()>1)
1728 ++fHeader->fNSelPrimTr1;
1729 if(track->Pt()>2)
1730 ++fHeader->fNSelPrimTr2;
1731 }
5fe1ca23 1732
1733 fHeader->fNCells = 0;
c0d2e1f2 1734 fHeader->fNCells0 = 0;
1735 fHeader->fNCells01 = 0;
1736 fHeader->fNCells03 = 0;
5fe1ca23 1737 fHeader->fNCells1 = 0;
1738 fHeader->fNCells2 = 0;
1739 fHeader->fNCells5 = 0;
1740 fHeader->fMaxCellE = 0;
1741
1742 AliVCaloCells *cells = fEsdCells;
1743 if (!cells)
1744 cells = fAodCells;
1745
1746 if (cells) {
1747 Int_t ncells = cells->GetNumberOfCells();
1748 for(Int_t j=0; j<ncells; ++j) {
1749 Double_t cellen = cells->GetAmplitude(j);
c0d2e1f2 1750 if (cellen>0.045)
1751 ++fHeader->fNCells0;
1752 if (cellen>0.1)
1753 ++fHeader->fNCells01;
1754 if (cellen>0.3)
1755 ++fHeader->fNCells03;
5fe1ca23 1756 if (cellen>1)
1757 ++fHeader->fNCells1;
1758 if (cellen>2)
1759 ++fHeader->fNCells2;
1760 if (cellen>5)
1761 ++fHeader->fNCells5;
1762 if (cellen>fHeader->fMaxCellE)
1763 fHeader->fMaxCellE = cellen;
1764 }
1765 fHeader->fNCells = ncells;
1766 }
1767
1768 fHeader->fNClus = 0;
1769 fHeader->fNClus1 = 0;
1770 fHeader->fNClus2 = 0;
1771 fHeader->fNClus5 = 0;
1772 fHeader->fMaxClusE = 0;
1773
1774 TObjArray *clusters = fEsdClusters;
1775 if (!clusters)
1776 clusters = fAodClusters;
1777
1778 if (clusters) {
1779 Int_t nclus = clusters->GetEntries();
1780 for(Int_t j=0; j<nclus; ++j) {
1781 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
6d4faab0 1782 if (!clus->IsEMCAL())
1783 continue;
5fe1ca23 1784 Double_t clusen = clus->E();
1785 if (clusen>1)
1786 ++fHeader->fNClus1;
1787 if (clusen>2)
1788 ++fHeader->fNClus2;
1789 if (clusen>5)
1790 ++fHeader->fNClus5;
1791 if (clusen>fHeader->fMaxClusE)
1792 fHeader->fMaxClusE = clusen;
1793 }
1794 fHeader->fNClus = nclus;
1795 }
1796
788ca675 1797 if (fAodEv) {
1798 am->LoadBranch("vertices");
1799 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1800 FillVertex(fPrimVert, pv);
1801 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1802 FillVertex(fSpdVert, sv);
1803 } else {
1804 am->LoadBranch("PrimaryVertex.");
1805 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1806 FillVertex(fPrimVert, pv);
1807 am->LoadBranch("SPDVertex.");
1808 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1809 FillVertex(fSpdVert, sv);
1810 am->LoadBranch("TPCVertex.");
1811 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1812 FillVertex(fTpcVert, tv);
fa443410 1813 }
788ca675 1814
788ca675 1815 fNtuple->Fill();
286b47a5 1816}
1817
1818//________________________________________________________________________
1819void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1820{
1821 // Fill histograms related to pions.
286b47a5 1822
296ea9b4 1823 Double_t vertex[3] = {0};
fa443410 1824 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1825
90d5b88b 1826 TObjArray *clusters = fEsdClusters;
1827 if (!clusters)
1828 clusters = fAodClusters;
fa443410 1829
90d5b88b 1830 if (clusters) {
1831 TLorentzVector clusterVec1;
1832 TLorentzVector clusterVec2;
1833 TLorentzVector pionVec;
1834
1835 Int_t nclus = clusters->GetEntries();
fa443410 1836 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1837 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1838 if (!clus1)
1839 continue;
1840 if (!clus1->IsEMCAL())
1841 continue;
3c661da5 1842 if (clus1->E()<fMinE)
6eb6260e 1843 continue;
a49742b5 1844 if (clus1->GetNCells()<fNminCells)
1845 continue;
f224d35b 1846 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1847 continue;
3c661da5 1848 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1849 continue;
fa443410 1850 clus1->GetMomentum(clusterVec1,vertex);
1851 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1852 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1853 if (!clus2)
1854 continue;
1855 if (!clus2->IsEMCAL())
1856 continue;
3c661da5 1857 if (clus2->E()<fMinE)
6eb6260e 1858 continue;
a49742b5 1859 if (clus2->GetNCells()<fNminCells)
1860 continue;
f224d35b 1861 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1862 continue;
3c661da5 1863 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1864 continue;
fa443410 1865 clus2->GetMomentum(clusterVec2,vertex);
1866 pionVec = clusterVec1 + clusterVec2;
1867 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1868 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1869 if (pionZgg < fAsymMax) {
1870 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1871 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1872 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1873 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1874 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1875 fHPionInvMasses[bin]->Fill(pionVec.M());
1876 }
fa443410 1877 }
1878 }
90d5b88b 1879 }
fa443410 1880}
d595acbb 1881
38727e64 1882//________________________________________________________________________
1883void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1884{
1885 // Fill additional MC information histograms.
1886
cfd7d5b2 1887 if (!fMcParts)
38727e64 1888 return;
1889
1890 // check if aod or esd mc mode and the fill histos
1891}
1892
6eb6260e 1893//________________________________________________________________________
323834f0 1894void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1895{
788ca675 1896 // Fill other histograms.
b5ba4932 1897
1898 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1899 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1900 if(!track)
1901 continue;
1902 fHPrimTrackPt->Fill(track->Pt());
1903 fHPrimTrackEta->Fill(track->Eta());
1904 fHPrimTrackPhi->Fill(track->Phi());
1905 }
788ca675 1906}
1907
f5e0f1e2 1908//________________________________________________________________________
1909void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1910{
1911 // Fill track histograms.
1912
1913 if (fSelPrimTracks) {
1914 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1915 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1916 if(!track)
1917 continue;
1918 fHPrimTrackPt->Fill(track->Pt());
1919 fHPrimTrackEta->Fill(track->Eta());
1920 fHPrimTrackPhi->Fill(track->Phi());
1921 }
1922 }
1923}
1924
788ca675 1925//__________________________________________________________________________________________________
1926void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1927{
1928 // Fill vertex from ESD vertex info.
1929
1930 v->fVx = esdv->GetX();
1931 v->fVy = esdv->GetY();
1932 v->fVz = esdv->GetZ();
1933 v->fVc = esdv->GetNContributors();
1934 v->fDisp = esdv->GetDispersion();
1935 v->fZres = esdv->GetZRes();
1936 v->fChi2 = esdv->GetChi2();
1937 v->fSt = esdv->GetStatus();
1938 v->fIs3D = esdv->IsFromVertexer3D();
1939 v->fIsZ = esdv->IsFromVertexerZ();
1940}
1941
1942//__________________________________________________________________________________________________
1943void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1944{
1945 // Fill vertex from AOD vertex info.
1946
1947 v->fVx = aodv->GetX();
1948 v->fVy = aodv->GetY();
1949 v->fVz = aodv->GetZ();
1950 v->fVc = aodv->GetNContributors();
1951 v->fChi2 = aodv->GetChi2();
6eb6260e 1952}
1953
d595acbb 1954//________________________________________________________________________
296ea9b4 1955Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1956{
1957 // Compute isolation based on cell content.
1958
1959 AliVCaloCells *cells = fEsdCells;
1960 if (!cells)
1961 cells = fAodCells;
1962 if (!cells)
1963 return 0;
1964
1965 Double_t cellIsolation = 0;
1966 Double_t rad2 = radius*radius;
1967 Int_t ncells = cells->GetNumberOfCells();
1968 for (Int_t i = 0; i<ncells; ++i) {
1969 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1970 Float_t eta=-1, phi=-1;
296ea9b4 1971 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1972 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1973 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1974 if(dist>rad2)
1975 continue;
0fbe8d4f 1976 Double_t cellE = cells->GetAmplitude(i);
296ea9b4 1977 cellIsolation += cellE;
1978 }
1979 return cellIsolation;
1980}
1981
1982//________________________________________________________________________
0fbe8d4f 1983Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
1984{
1985 // Get maximum energy of attached cell.
1986
1987 Double_t ret = 0;
1988 Int_t ncells = cluster->GetNCells();
1989 if (fEsdCells) {
1990 for (Int_t i=0; i<ncells; i++) {
1991 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1992 ret += e;
1993 }
1994 } else {
1995 for (Int_t i=0; i<ncells; i++) {
1996 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1997 ret += e;
1998 }
1999 }
2000 return ret;
2001}
2002
2003//________________________________________________________________________
2004Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 2005{
90d5b88b 2006 // Get maximum energy of attached cell.
2007
788ca675 2008 id = -1;
90d5b88b 2009 Double_t maxe = 0;
90d5b88b 2010 Int_t ncells = cluster->GetNCells();
2011 if (fEsdCells) {
2012 for (Int_t i=0; i<ncells; i++) {
27422847 2013 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 2014 if (e>maxe) {
90d5b88b 2015 maxe = e;
788ca675 2016 id = cluster->GetCellAbsId(i);
f0897c18 2017 }
90d5b88b 2018 }
2019 } else {
2020 for (Int_t i=0; i<ncells; i++) {
27422847 2021 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 2022 if (e>maxe)
2023 maxe = e;
788ca675 2024 id = cluster->GetCellAbsId(i);
90d5b88b 2025 }
2026 }
6eb6260e 2027 return maxe;
d595acbb 2028}
2029
2030//________________________________________________________________________
0fbe8d4f 2031void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 2032{
90d5b88b 2033 // Calculate the (E) weighted variance along the longer (eigen) axis.
2034
6bf90832 2035 sigmaMax = 0; // cluster variance along its longer axis
2036 sigmaMin = 0; // cluster variance along its shorter axis
2037 Double_t Ec = c->E(); // cluster energy
296ea9b4 2038 if(Ec<=0)
2039 return;
6bf90832 2040 Double_t Xc = 0 ; // cluster first moment along X
2041 Double_t Yc = 0 ; // cluster first moment along Y
2042 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
2043 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
2044 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 2045
2046 AliVCaloCells *cells = fEsdCells;
2047 if (!cells)
2048 cells = fAodCells;
2049
6eb6260e 2050 if (!cells)
2051 return;
2052
6bf90832 2053 Int_t ncells = c->GetNCells();
6eb6260e 2054 if (ncells==1)
2055 return;
2056
6eb6260e 2057 for(Int_t j=0; j<ncells; ++j) {
6bf90832 2058 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 2059 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 2060 TVector3 pos;
2061 fGeom->GetGlobal(id,pos);
6eb6260e 2062 Xc += cellen*pos.X();
2063 Yc += cellen*pos.Y();
2064 Sxx += cellen*pos.X()*pos.X();
2065 Syy += cellen*pos.Y()*pos.Y();
2066 Sxy += cellen*pos.X()*pos.Y();
2067 }
2068 Xc /= Ec;
2069 Yc /= Ec;
2070 Sxx /= Ec;
2071 Syy /= Ec;
2072 Sxy /= Ec;
2073 Sxx -= Xc*Xc;
2074 Syy -= Yc*Yc;
2075 Sxy -= Xc*Yc;
f0897c18 2076 Sxx = TMath::Abs(Sxx);
2077 Syy = TMath::Abs(Syy);
296ea9b4 2078 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2079 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2080 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2081 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 2082}
90d5b88b 2083
6bf90832 2084//________________________________________________________________________
0fbe8d4f 2085Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 2086{
2087 // Calculate number of attached cells above emin.
2088
6bf90832 2089 AliVCaloCells *cells = fEsdCells;
2090 if (!cells)
2091 cells = fAodCells;
6bf90832 2092 if (!cells)
296ea9b4 2093 return 0;
6bf90832 2094
296ea9b4 2095 Int_t n = 0;
6bf90832 2096 Int_t ncells = c->GetNCells();
2097 for(Int_t j=0; j<ncells; ++j) {
2098 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2099 Double_t cellen = cells->GetCellAmplitude(id);
2100 if (cellen>=emin)
2101 ++n;
2102 }
2103 return n;
2104}
2105
296ea9b4 2106//________________________________________________________________________
b6c599fe 2107Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 2108{
2109 // Compute isolation based on tracks.
2110
2111 Double_t trkIsolation = 0;
2112 Double_t rad2 = radius*radius;
b6c599fe 2113 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 2114 for(Int_t j = 0; j<ntrks; ++j) {
2115 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
2116 if (!track)
2117 continue;
b6c599fe 2118 if (track->Pt()<pt)
2119 continue;
296ea9b4 2120 Float_t eta = track->Eta();
2121 Float_t phi = track->Phi();
0fbe8d4f 2122 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 2123 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2124 if(dist>rad2)
2125 continue;
2126 trkIsolation += track->Pt();
2127 }
2128 return trkIsolation;
2129}
3a952328 2130
0fbe8d4f 2131//________________________________________________________________________
2132Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2133{
2134 // Returns if cluster shared across super modules.
2135
2136 AliVCaloCells *cells = fEsdCells;
2137 if (!cells)
2138 cells = fAodCells;
2139 if (!cells)
2140 return 0;
2141
2142 Int_t n = -1;
2143 Int_t ncells = c->GetNCells();
2144 for(Int_t j=0; j<ncells; ++j) {
2145 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2146 Int_t got = id / (24*48);
2147 if (n==-1) {
2148 n = got;
2149 continue;
2150 }
2151 if (got!=n)
2152 return 1;
2153 }
2154 return 0;
2155}
3a952328 2156
38727e64 2157//________________________________________________________________________
807016ea 2158void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
38727e64 2159{
56fd6cb2 2160 // Print recursively daughter information.
2161
2162 if (!p || !arr)
2163 return;
2164
807016ea 2165 const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2166 if (!amc)
2167 return;
2168 for (Int_t i=0; i<level; ++i) printf(" ");
2169 amc->Print();
2170
2171 Int_t n = amc->GetNDaughters();
2172 for (Int_t i=0; i<n; ++i) {
2173 Int_t d = amc->GetDaughter(i);
2174 const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2175 PrintDaughters(dmc,arr,level+1);
2176 }
2177}
2178
2179//________________________________________________________________________
2180void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2181{
2182 // Print recursively daughter information.
2183
2184 if (!p || !arr)
2185 return;
2186
2187 for (Int_t i=0; i<level; ++i) printf(" ");
2188 Int_t d1 = p->GetFirstDaughter();
2189 Int_t d2 = p->GetLastDaughter();
2190 printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2191 p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2192 if (d1<0)
2193 return;
2194 if (d2<0)
2195 d2=d1;
2196 for (Int_t i=d1;i<=d2;++i) {
2197 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
6d4faab0 2198 PrintDaughters(dmc,arr,level+1);
56fd6cb2 2199 }
38727e64 2200}
2201
2202//________________________________________________________________________
2203void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2204{
56fd6cb2 2205 // Print track ref array.
2206
2207 if (!p)
2208 return;
2209
2210 Int_t n = p->GetNumberOfTrackReferences();
2211 for (Int_t i=0; i<n; ++i) {
2212 AliTrackReference *ref = p->GetTrackReference(i);
bc8a45bd 2213 if (!ref)
2214 continue;
56fd6cb2 2215 ref->SetUserId(ref->DetectorId());
2216 ref->Print();
2217 }
38727e64 2218}
2219
807016ea 2220//________________________________________________________________________
2221void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2222{
2223 // Process and create daughters.
2224
2225 if (!p || !arr)
2226 return;
2227
2228 AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2229 if (!amc)
2230 return;
2231
2232 //amc->Print();
2233
2234 Int_t nparts = arr->GetEntries();
2235 Int_t nents = fMcParts->GetEntries();
2236
2237 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2238 newp->fPt = amc->Pt();
2239 newp->fEta = amc->Eta();
2240 newp->fPhi = amc->Phi();
2241 if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2242 TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2243 newp->fVR = vec.Perp();
2244 newp->fVEta = vec.Eta();
2245 newp->fVPhi = vec.Phi();
2246 }
2247 newp->fPid = amc->PdgCode();
8c56d760 2248 newp->fLab = nents;
807016ea 2249 Int_t moi = amc->GetMother();
2250 if (moi>=0&&moi<nparts) {
2251 const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2252 moi = mmc->GetUniqueID();
2253 }
2254 newp->fMo = moi;
2255 p->SetUniqueID(nents);
8c56d760 2256
2257 // TODO: Determine which detector was hit
2258 //newp->fDet = ???
2259
807016ea 2260 Int_t n = amc->GetNDaughters();
2261 for (Int_t i=0; i<n; ++i) {
2262 Int_t d = amc->GetDaughter(i);
2263 if (d<=index || d>=nparts)
2264 continue;
2265 AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2266 ProcessDaughters(dmc,d,arr);
2267 }
2268}
2269
2270//________________________________________________________________________
2271void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2272{
2273 // Process and create daughters.
2274
2275 if (!p || !arr)
2276 return;
2277
2278 Int_t d1 = p->GetFirstDaughter();
2279 Int_t d2 = p->GetLastDaughter();
2280 if (0) {
2281 printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2282 index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
bc8a45bd 2283 PrintTrackRefs(p);
807016ea 2284 }
2285 Int_t nents = fMcParts->GetEntries();
2286
2287 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2288 newp->fPt = p->Pt();
2289 newp->fEta = p->Eta();
2290 newp->fPhi = p->Phi();
2291 if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2292 TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2293 newp->fVR = vec.Perp();
2294 newp->fVEta = vec.Eta();
2295 newp->fVPhi = vec.Phi();
2296 }
2297 newp->fPid = p->PdgCode();
8c56d760 2298 newp->fLab = nents;
807016ea 2299 Int_t moi = p->GetMother();
2300 if (moi>=0) {
2301 const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2302 moi = mmc->GetUniqueID();
2303 }
2304 newp->fMo = moi;
2305 p->SetUniqueID(nents);
2306
2307 Int_t nref = p->GetNumberOfTrackReferences();
2308 if (nref>0) {
2309 AliTrackReference *ref = p->GetTrackReference(nref-1);
2310 if (ref) {
2311 newp->fDet = ref->DetectorId();
2312 }
2313 }
2314
2315 if (d1<0)
2316 return;
2317 if (d2<0)
2318 d2=d1;
2319 for (Int_t i=d1;i<=d2;++i) {
2320 AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2321 if (dmc->P()<0.01)
2322 continue;
2323 ProcessDaughters(dmc,i,arr);
2324 }
2325}