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