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