]> 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
Added more severe checks if emcal matrizes are there
[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"
18#include "AliAODVertex.h"
19#include "AliAnalysisManager.h"
20#include "AliAnalysisTaskEMCALClusterizeFast.h"
296ea9b4 21#include "AliCDBManager.h"
ea3fd2d5 22#include "AliCentrality.h"
ea3fd2d5 23#include "AliEMCALGeoUtils.h"
0ec74551 24#include "AliEMCALGeometry.h"
0fbe8d4f 25#include "AliEMCALRecPoint.h"
296ea9b4 26#include "AliEMCALRecoUtils.h"
3a952328 27#include "AliESDCaloTrigger.h"
717fe7de 28#include "AliESDEvent.h"
5fe1ca23 29#include "AliESDUtils.h"
717fe7de 30#include "AliESDVertex.h"
296ea9b4 31#include "AliESDtrack.h"
0ec74551 32#include "AliESDtrackCuts.h"
b6c599fe 33#include "AliEventplane.h"
296ea9b4 34#include "AliGeomManager.h"
b3ee6797 35#include "AliInputEventHandler.h"
43cfaa06 36#include "AliLog.h"
296ea9b4 37#include "AliMagF.h"
5fe1ca23 38#include "AliMultiplicity.h"
0ec74551 39#include "AliTrackerBase.h"
ea3fd2d5 40
41ClassImp(AliAnalysisTaskEMCALPi0PbPb)
717fe7de 42
ea3fd2d5 43//________________________________________________________________________
44AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
45 : AliAnalysisTaskSE(name),
717fe7de 46 fCentVar("V0M"),
47 fCentFrom(0),
48 fCentTo(100),
76332037 49 fVtxZMin(-10),
50 fVtxZMax(+10),
43cfaa06 51 fUseQualFlag(1),
717fe7de 52 fClusName(),
f5d4ab70 53 fDoNtuple(0),
a49742b5 54 fDoAfterburner(0),
f224d35b 55 fAsymMax(1),
a49742b5 56 fNminCells(1),
3c661da5 57 fMinE(0.100),
f224d35b 58 fMinErat(0),
59 fMinEcc(0),
6bf90832 60 fGeoName("EMCAL_FIRSTYEARV1"),
b6c599fe 61 fMinNClusPerTr(50),
296ea9b4 62 fIsoDist(0.2),
b3ee6797 63 fTrClassNames(""),
0ec74551 64 fTrCuts(0),
b6c599fe 65 fPrimTrCuts(0),
66 fDoTrMatGeom(0),
3a952328 67 fTrainMode(0),
68 fMarkCells(),
69 fMinL0Time(-1),
70 fMaxL0Time(1024),
27c2e3d9 71 fIsGeoMatsSet(0),
d9f26424 72 fNEvs(0),
d595acbb 73 fGeom(0),
296ea9b4 74 fReco(0),
717fe7de 75 fOutput(0),
b3ee6797 76 fTrClassNamesArr(0),
717fe7de 77 fEsdEv(0),
78 fAodEv(0),
43cfaa06 79 fRecPoints(0),
717fe7de 80 fEsdClusters(0),
81 fEsdCells(0),
82 fAodClusters(0),
286b47a5 83 fAodCells(0),
fa443410 84 fPtRanges(0),
296ea9b4 85 fSelTracks(0),
b6c599fe 86 fSelPrimTracks(0),
3a952328 87 fNAmpInTrigger(0),
88 fAmpInTrigger(0),
788ca675 89 fNtuple(0),
90 fHeader(0),
91 fPrimVert(0),
92 fSpdVert(0),
93 fTpcVert(0),
94 fClusters(0),
3a952328 95 fTriggers(0),
fa443410 96 fHCuts(0x0),
97 fHVertexZ(0x0),
76332037 98 fHVertexZ2(0x0),
fa443410 99 fHCent(0x0),
100 fHCentQual(0x0),
b3ee6797 101 fHTclsBeforeCuts(0x0),
102 fHTclsAfterCuts(0x0),
d595acbb 103 fHColuRow(0x0),
104 fHColuRowE(0x0),
105 fHCellMult(0x0),
106 fHCellE(0x0),
107 fHCellH(0x0),
6eb6260e 108 fHCellM(0x0),
a49742b5 109 fHCellM2(0x0),
296ea9b4 110 fHCellFreqNoCut(0x0),
2e4d8148 111 fHCellFreqCut100M(0x0),
112 fHCellFreqCut300M(0x0),
113 fHCellFreqE(0x0),
296ea9b4 114 fHCellCheckE(0x0),
6eb6260e 115 fHClustEccentricity(0),
fa443410 116 fHClustEtaPhi(0x0),
117 fHClustEnergyPt(0x0),
118 fHClustEnergySigma(0x0),
d595acbb 119 fHClustSigmaSigma(0x0),
6eb6260e 120 fHClustNCellEnergyRatio(0x0),
b6c599fe 121 fHMatchDr(0x0),
122 fHMatchDz(0x0),
123 fHMatchEp(0x0),
fa443410 124 fHPionEtaPhi(0x0),
125 fHPionMggPt(0x0),
6eb6260e 126 fHPionMggAsym(0x0),
127 fHPionMggDgg(0x0)
ea3fd2d5 128{
717fe7de 129 // Constructor.
ea3fd2d5 130
d595acbb 131 if (!name)
132 return;
133 SetName(name);
ea3fd2d5 134 DefineInput(0, TChain::Class());
ea3fd2d5 135 DefineOutput(1, TList::Class());
3a952328 136 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks EMCALTrigger."
296ea9b4 137 "AOD:header,vertices,emcalCells,tracks";
ea3fd2d5 138}
ea3fd2d5 139
717fe7de 140//________________________________________________________________________
141AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
142{
143 // Destructor.
ea3fd2d5 144
b3ee6797 145 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
146 delete fOutput; fOutput = 0;
147 }
fa443410 148 delete fPtRanges; fPtRanges = 0;
d595acbb 149 delete fGeom; fGeom = 0;
296ea9b4 150 delete fReco; fReco = 0;
b3ee6797 151 delete fTrClassNamesArr;
296ea9b4 152 delete fSelTracks;
b6c599fe 153 delete fSelPrimTracks;
3a952328 154 delete [] fAmpInTrigger;
d595acbb 155 delete [] fHColuRow;
156 delete [] fHColuRowE;
157 delete [] fHCellMult;
296ea9b4 158 delete [] fHCellFreqNoCut;
2e4d8148 159 delete [] fHCellFreqCut100M;
160 delete [] fHCellFreqCut300M;
ea3fd2d5 161}
717fe7de 162
ea3fd2d5 163//________________________________________________________________________
164void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
165{
717fe7de 166 // Create user objects here.
ea3fd2d5 167
f3582e89 168 cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
169 cout << " fCentVar: " << fCentVar << endl;
170 cout << " fCentFrom: " << fCentFrom << endl;
171 cout << " fCentTo: " << fCentTo << endl;
172 cout << " fVtxZMin: " << fVtxZMin << endl;
173 cout << " fVtxZMax: " << fVtxZMax << endl;
174 cout << " fUseQualFlag: " << fUseQualFlag << endl;
175 cout << " fClusName: \"" << fClusName << "\"" << endl;
176 cout << " fDoNtuple: " << fDoNtuple << endl;
177 cout << " fDoAfterburner: " << fDoAfterburner << endl;
178 cout << " fAsymMax: " << fAsymMax << endl;
179 cout << " fNminCells: " << fNminCells << endl;
180 cout << " fMinE: " << fMinE << endl;
181 cout << " fMinErat: " << fMinErat << endl;
182 cout << " fMinEcc: " << fMinEcc << endl;
183 cout << " fGeoName: \"" << fGeoName << "\"" << endl;
184 cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
185 cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
f3582e89 186 cout << " fTrCuts: " << fTrCuts << endl;
187 cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
188 cout << " fDoTrMatGeom: " << fDoTrMatGeom << endl;
189 cout << " fTrainMode: " << fTrainMode << endl;
190 cout << " fMarkCells: " << fMarkCells << endl;
191 cout << " fMinL0Time: " << fMinL0Time << endl;
192 cout << " fMaxL0Time: " << fMaxL0Time << endl;
193
296ea9b4 194 fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
195 fReco = new AliEMCALRecoUtils();
b3ee6797 196 fTrClassNamesArr = fTrClassNames.Tokenize(" ");
717fe7de 197 fOutput = new TList();
198 fOutput->SetOwner();
296ea9b4 199 fSelTracks = new TObjArray;
b6c599fe 200 fSelPrimTracks = new TObjArray;
ea3fd2d5 201
f5d4ab70 202 if (fDoNtuple) {
203 TFile *f = OpenFile(1);
6bf90832 204 if (f) {
205 f->SetCompressionLevel(2);
788ca675 206 fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
6bf90832 207 fNtuple->SetDirectory(f);
3a952328 208 if (fTrainMode) {
209 fNtuple->SetAutoFlush(-2*1024*1024);
210 fNtuple->SetAutoSave(0);
211 } else {
212 fNtuple->SetAutoFlush(-32*1024*1024);
213 fNtuple->SetAutoSave(0);
214 }
215
788ca675 216 fHeader = new AliStaHeader;
217 fNtuple->Branch("header", &fHeader, 16*1024, 99);
218 fPrimVert = new AliStaVertex;
219 fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
220 fSpdVert = new AliStaVertex;
221 fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
222 fTpcVert = new AliStaVertex;
223 fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
224 if (TClass::GetClass("AliStaCluster"))
225 TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
3a952328 226 fClusters = new TClonesArray("AliStaCluster");
788ca675 227 fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
3a952328 228 if (TClass::GetClass("AliStaTrigger"))
229 TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
230 fTriggers = new TClonesArray("AliStaTrigger");
231 fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
6bf90832 232 }
f5d4ab70 233 }
234
002dcebe 235 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 236 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
237 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
002dcebe 238
d595acbb 239 // histograms
a49742b5 240 TH1::SetDefaultSumw2(kTRUE);
241 TH2::SetDefaultSumw2(kTRUE);
b3ee6797 242 fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
243 fHCuts->GetXaxis()->SetBinLabel(1,"All");
244 fHCuts->GetXaxis()->SetBinLabel(2,"PS");
245 fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
246 fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
247 fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
fa443410 248 fOutput->Add(fHCuts);
249 fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
250 fHVertexZ->SetXTitle("z [cm]");
251 fOutput->Add(fHVertexZ);
76332037 252 fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
253 fHVertexZ2->SetXTitle("z [cm]");
254 fOutput->Add(fHVertexZ2);
f2b8fca6 255 fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
fa443410 256 fHCent->SetXTitle(fCentVar.Data());
76332037 257 fOutput->Add(fHCent);
f2b8fca6 258 fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
fa443410 259 fHCentQual->SetXTitle(fCentVar.Data());
260 fOutput->Add(fHCentQual);
2e4d8148 261 fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
262 fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
b3ee6797 263 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
264 const char *name = fTrClassNamesArr->At(i)->GetName();
265 fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
266 fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
267 }
268 fOutput->Add(fHTclsBeforeCuts);
269 fOutput->Add(fHTclsAfterCuts);
90d5b88b 270
d595acbb 271 // histograms for cells
296ea9b4 272 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
273 fHColuRow = new TH2*[nsm];
274 fHColuRowE = new TH2*[nsm];
275 fHCellMult = new TH1*[nsm];
276 for (Int_t i = 0; i<nsm; ++i) {
2e4d8148 277 fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
d595acbb 278 fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
279 fHColuRow[i]->SetXTitle("col (i#eta)");
280 fHColuRow[i]->SetYTitle("row (i#phi)");
2e4d8148 281 fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
90d5b88b 282 fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
d595acbb 283 fHColuRowE[i]->SetXTitle("col (i#eta)");
284 fHColuRowE[i]->SetYTitle("row (i#phi)");
285 fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
286 fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
90d5b88b 287 fHCellMult[i]->SetXTitle("# of cells");
d595acbb 288 fOutput->Add(fHColuRow[i]);
289 fOutput->Add(fHColuRowE[i]);
290 fOutput->Add(fHCellMult[i]);
291 }
2e4d8148 292 fHCellE = new TH1F("hCellE","",250,0.,25.);
d595acbb 293 fHCellE->SetXTitle("E_{cell} [GeV]");
294 fOutput->Add(fHCellE);
fb5a0b69 295 fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
6eb6260e 296 fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
d595acbb 297 fOutput->Add(fHCellH);
fb5a0b69 298 fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
6eb6260e 299 fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
300 fOutput->Add(fHCellM);
fb5a0b69 301 fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
f4ec884e 302 fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
a49742b5 303 fOutput->Add(fHCellM2);
90d5b88b 304
2e4d8148 305 fHCellFreqNoCut = new TH1*[nsm];
306 fHCellFreqCut100M = new TH1*[nsm];
307 fHCellFreqCut300M = new TH1*[nsm];
308 fHCellFreqE = new TH1*[nsm];
296ea9b4 309 for (Int_t i = 0; i<nsm; ++i){
310 Double_t lbin = i*24*48-0.5;
311 Double_t hbin = lbin+24*48;
2e4d8148 312 fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i),
313 Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin);
314 fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
315 Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
316 fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
317 Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
318 fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i),
319 Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
296ea9b4 320 fOutput->Add(fHCellFreqNoCut[i]);
2e4d8148 321 fOutput->Add(fHCellFreqCut100M[i]);
322 fOutput->Add(fHCellFreqCut300M[i]);
323 fOutput->Add(fHCellFreqE[i]);
296ea9b4 324 }
3a952328 325 if (!fMarkCells.IsNull()) {
296ea9b4 326 fHCellCheckE = new TH1*[24*48*nsm];
408dc04c 327 memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
3a952328 328 TObjArray *cells = fMarkCells.Tokenize(" ");
329 Int_t n = cells->GetEntries();
330 Int_t *tcs = new Int_t[n];
331 for (Int_t i=0;i<n;++i) {
332 TString name(cells->At(i)->GetName());
333 tcs[i]=name.Atoi();
334 }
335 for (Int_t i = 0; i<n; ++i) {
296ea9b4 336 Int_t c=tcs[i];
337 if (c<24*48*nsm) {
3a952328 338 fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
296ea9b4 339 fOutput->Add(fHCellCheckE[i]);
340 }
341 }
3a952328 342 delete cells;
343 delete [] tcs;
296ea9b4 344 }
345
d595acbb 346 // histograms for clusters
3a952328 347 if (!fTrainMode) {
348 fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
349 fHClustEccentricity->SetXTitle("#epsilon_{C}");
350 fOutput->Add(fHClustEccentricity);
351 fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
352 fHClustEtaPhi->SetXTitle("#eta");
353 fHClustEtaPhi->SetYTitle("#varphi");
354 fOutput->Add(fHClustEtaPhi);
355 fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
356 fHClustEnergyPt->SetXTitle("E [GeV]");
357 fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
358 fOutput->Add(fHClustEnergyPt);
359 fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
360 fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
361 fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
362 fOutput->Add(fHClustEnergySigma);
363 fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
364 fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
365 fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
366 fOutput->Add(fHClustSigmaSigma);
367 fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
368 fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
369 fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
370 fOutput->Add(fHClustNCellEnergyRatio);
371 }
90d5b88b 372
b6c599fe 373 // histograms for track matching
374 fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
375 fOutput->Add(fHMatchDr);
376 fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
377 fOutput->Add(fHMatchDz);
378 fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
379 fOutput->Add(fHMatchEp);
3a952328 380
d595acbb 381 // histograms for pion candidates
3a952328 382 if (!fTrainMode) {
383 fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
384 fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
385 fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
386 fOutput->Add(fHPionEtaPhi);
387 fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
388 fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
389 fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
390 fOutput->Add(fHPionMggPt);
391 fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
392 fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
393 fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
394 fOutput->Add(fHPionMggAsym);
395 fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
396 fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
397 fHPionMggDgg->SetYTitle("opening angle [grad]");
398 fOutput->Add(fHPionMggDgg);
399 const Int_t nbins = 20;
400 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};
401 fPtRanges = new TAxis(nbins-1,xbins);
402 for (Int_t i = 0; i<=nbins; ++i) {
403 fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
404 fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
405 if (i==0)
406 fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
407 else if (i==nbins)
408 fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
409 else
410 fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
411 fOutput->Add(fHPionInvMasses[i]);
412 }
fa443410 413 }
296ea9b4 414
717fe7de 415 PostData(1, fOutput);
ea3fd2d5 416}
417
418//________________________________________________________________________
419void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
420{
717fe7de 421 // Called for each event.
422
423 if (!InputEvent())
424 return;
425
43cfaa06 426 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
427 fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
44310bce 428 UInt_t offtrigger = 0;
43cfaa06 429 if (fEsdEv) {
430 am->LoadBranch("AliESDRun.");
431 am->LoadBranch("AliESDHeader.");
9809ed8c 432 UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
433 UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
434 Bool_t desc1 = (mask1 >> 18) & 0x1;
435 Bool_t desc2 = (mask2 >> 18) & 0x1;
436 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
437 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
438 mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
439 mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
44310bce 440 return;
441 }
442 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
43cfaa06 443 } else {
444 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
408dc04c 445 if (!fAodEv) {
446 AliFatal("Neither ESD nor AOD event found");
447 return;
448 }
43cfaa06 449 am->LoadBranch("header");
44310bce 450 offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
451 }
452 if (offtrigger & AliVEvent::kFastOnly) {
453 AliWarning(Form("EMCAL not in fast only partition"));
454 return;
43cfaa06 455 }
b6c599fe 456 if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
27c2e3d9 457 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
458 AliCDBManager *cdb = AliCDBManager::Instance();
459 if (!cdb->IsDefaultStorageSet())
460 cdb->SetDefaultStorage("raw://");
461 Int_t runno = InputEvent()->GetRunNumber();
462 if (runno != cdb->GetRun())
463 cdb->SetRun(runno);
464 AliGeomManager::LoadGeometry();
465 }
466
467 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
468 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
9809ed8c 469 for (Int_t i=0; i<nsm; ++i) {
470 const TGeoHMatrix *geom = 0;
471 if (fEsdEv)
472 geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
473 else
474 geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
475 if (!geom)
476 continue;
477 geom->Print();
478 fGeom->SetMisalMatrix(geom,i);
27c2e3d9 479 }
480 fIsGeoMatsSet = kTRUE;
481 }
482
483 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
484 if (fEsdEv) {
485 const AliESDRun *erun = fEsdEv->GetESDRun();
486 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
487 erun->GetCurrentDip(),
488 AliMagF::kConvLHC,
489 kFALSE,
490 erun->GetBeamEnergy(),
491 erun->GetBeamType());
492 TGeoGlobalMagField::Instance()->SetField(field);
493 } else {
494 Double_t pol = -1; //polarity
495 Double_t be = -1; //beam energy
496 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
497 Int_t runno = fAodEv->GetRunNumber();
498 if (runno>=136851 && runno<138275) {
499 pol = -1;
500 be = 2760;
501 btype = AliMagF::kBeamTypeAA;
502 } else if (runno>=138275 && runno<=139517) {
503 pol = +1;
504 be = 2760;
505 btype = AliMagF::kBeamTypeAA;
506 } else {
507 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
508 }
509 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
510 }
511 }
512
b3ee6797 513 Int_t cut = 1;
514 fHCuts->Fill(cut++);
515
516 TString trgclasses;
517 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
518 if (h) {
519 trgclasses = fEsdEv->GetFiredTriggerClasses();
520 } else {
521 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
522 if (h2)
523 trgclasses = h2->GetFiredTriggerClasses();
524 }
525 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
526 const char *name = fTrClassNamesArr->At(i)->GetName();
527 if (trgclasses.Contains(name))
528 fHTclsBeforeCuts->Fill(1+i);
529 }
530
531 UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
532 if (res==0)
533 return;
534
fa443410 535 fHCuts->Fill(cut++);
286b47a5 536
717fe7de 537 const AliCentrality *centP = InputEvent()->GetCentrality();
538 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
fa443410 539 fHCent->Fill(cent);
717fe7de 540 if (cent<fCentFrom||cent>fCentTo)
286b47a5 541 return;
43cfaa06 542
fa443410 543 fHCuts->Fill(cut++);
286b47a5 544
43cfaa06 545 if (fUseQualFlag) {
546 if (centP->GetQuality()>0)
547 return;
548 }
286b47a5 549
fa443410 550 fHCentQual->Fill(cent);
551 fHCuts->Fill(cut++);
717fe7de 552
43cfaa06 553 if (fEsdEv) {
fa443410 554 am->LoadBranch("PrimaryVertex.");
555 am->LoadBranch("SPDVertex.");
556 am->LoadBranch("TPCVertex.");
43cfaa06 557 } else {
558 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
559 am->LoadBranch("vertices");
3c661da5 560 if (!fAodEv) return;
43cfaa06 561 }
562
563 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
717fe7de 564 if (!vertex)
565 return;
566
fa443410 567 fHVertexZ->Fill(vertex->GetZ());
286b47a5 568
717fe7de 569 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
570 return;
286b47a5 571
fa443410 572 fHCuts->Fill(cut++);
76332037 573 fHVertexZ2->Fill(vertex->GetZ());
717fe7de 574
b3ee6797 575 // count number of accepted events
576 ++fNEvs;
577
578 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
579 const char *name = fTrClassNamesArr->At(i)->GetName();
580 if (trgclasses.Contains(name))
581 fHTclsAfterCuts->Fill(1+i);
582 }
583
43cfaa06 584 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
b6c599fe 585 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 586 fEsdCells = 0; // will be set if ESD input used
b6c599fe 587 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 588 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
589 fAodCells = 0; // will be set if AOD input used
717fe7de 590
43cfaa06 591 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
592 Bool_t clusattached = 0;
593 Bool_t recalibrated = 0;
594 if (1 && !fClusName.IsNull()) {
595 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
596 TObjArray *ts = am->GetTasks();
597 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
598 if (cltask && cltask->GetClusters()) {
d595acbb 599 fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
43cfaa06 600 clusattached = cltask->GetAttachClusters();
601 if (cltask->GetCalibData()!=0)
602 recalibrated = kTRUE;
603 }
604 }
605 if (1 && AODEvent() && !fClusName.IsNull()) {
717fe7de 606 TList *l = AODEvent()->GetList();
607 TClonesArray *clus = 0;
608 if (l) {
609 clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
610 fAodClusters = clus;
611 }
ea3fd2d5 612 }
717fe7de 613
614 if (fEsdEv) { // ESD input mode
43cfaa06 615 if (1 && (!fRecPoints||clusattached)) {
616 if (!clusattached)
617 am->LoadBranch("CaloClusters");
717fe7de 618 TList *l = fEsdEv->GetList();
619 TClonesArray *clus = 0;
620 if (l) {
621 clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
622 fEsdClusters = clus;
ea3fd2d5 623 }
624 }
43cfaa06 625 if (1) {
626 if (!recalibrated)
627 am->LoadBranch("EMCALCells.");
717fe7de 628 fEsdCells = fEsdEv->GetEMCALCells();
629 }
630 } else if (fAodEv) { // AOD input mode
43cfaa06 631 if (1 && (!fAodClusters || clusattached)) {
632 if (!clusattached)
633 am->LoadBranch("caloClusters");
717fe7de 634 TList *l = fAodEv->GetList();
635 TClonesArray *clus = 0;
636 if (l) {
637 clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
638 fAodClusters = clus;
ea3fd2d5 639 }
640 }
43cfaa06 641 if (1) {
642 if (!recalibrated)
643 am->LoadBranch("emcalCells");
717fe7de 644 fAodCells = fAodEv->GetEMCALCells();
645 }
646 } else {
647 AliFatal("Impossible to not have either pointer to ESD or AOD event");
648 }
ea3fd2d5 649
43cfaa06 650 if (1) {
651 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
652 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
653 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
654 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
655 AliDebug(2,Form("fAodCells set: %p", fAodCells));
656 }
657
a49742b5 658 if (fDoAfterburner)
659 ClusterAfterburner();
6eb6260e 660
3a952328 661 CalcCaloTriggers();
b6c599fe 662 CalcPrimTracks();
3a952328 663 CalcTracks();
296ea9b4 664 CalcClusterProps();
665
286b47a5 666 FillCellHists();
3a952328 667 if (!fTrainMode) {
668 FillClusHists();
669 FillPionHists();
670 FillOtherHists();
671 }
788ca675 672 FillNtuple();
ea3fd2d5 673
3a952328 674 if (fTrainMode) {
675 fSelTracks->Clear();
676 fSelPrimTracks->Clear();
677 fTriggers->Clear();
678 }
de4cee41 679
717fe7de 680 PostData(1, fOutput);
ea3fd2d5 681}
682
683//________________________________________________________________________
684void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
685{
717fe7de 686 // Terminate called at the end of analysis.
f5d4ab70 687
688 if (fNtuple) {
689 TFile *f = OpenFile(1);
690 if (f)
691 fNtuple->Write();
692 }
d9f26424 693
b3ee6797 694 AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
ea3fd2d5 695}
ea3fd2d5 696
296ea9b4 697//________________________________________________________________________
3a952328 698void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
296ea9b4 699{
3a952328 700 // Calculate triggers
296ea9b4 701
3a952328 702 fTriggers->Clear();
296ea9b4 703
3a952328 704 AliVCaloCells *cells = fEsdCells;
705 if (!cells)
706 cells = fAodCells;
707 if (!cells)
708 return;
709
710 Int_t ncells = cells->GetNumberOfCells();
711 if (ncells<=0)
712 return;
713
714 if (ncells>fNAmpInTrigger) {
715 delete [] fAmpInTrigger;
716 fAmpInTrigger = new Float_t[ncells];
717 fNAmpInTrigger = ncells;
296ea9b4 718 }
3a952328 719 for (Int_t i=0;i<ncells;++i)
720 fAmpInTrigger[i] = 0;
296ea9b4 721
3a952328 722 std::map<Short_t,Short_t> map;
723 for (Short_t pos=0;pos<ncells;++pos) {
724 Short_t id = cells->GetCellNumber(pos);
725 map[id]=pos;
726 }
727
728 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
729 am->LoadBranch("EMCALTrigger.");
730
731 AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
732 if (!triggers)
733 return;
734 if (triggers->GetEntries()<=0)
b6c599fe 735 return;
736
3a952328 737 triggers->Reset();
738 Int_t ntrigs=0;
739 while (triggers->Next()) {
740 Int_t gCol=0, gRow=0, ntimes=0;
741 triggers->GetPosition(gCol,gRow);
742 triggers->GetNL0Times(ntimes);
743 if (ntimes<1)
744 continue;
745 Float_t amp=0;
746 triggers->GetAmplitude(amp);
747 Int_t find = -1;
748 fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
749 if (find<0)
750 continue;
751 Int_t cidx[4] = {-1};
752 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
753 if (!ret)
754 continue;
755 Int_t trgtimes[25];
756 triggers->GetL0Times(trgtimes);
757 Int_t mintime = trgtimes[0];
758 Int_t maxtime = trgtimes[0];
759 Bool_t trigInTimeWindow = 0;
760 for (Int_t i=0;i<ntimes;++i) {
761 if (trgtimes[i]<mintime)
762 mintime = trgtimes[i];
763 if (maxtime<trgtimes[i])
764 maxtime = trgtimes[i];
f2b8fca6 765 if ((fMinL0Time<=trgtimes[i]) && (fMaxL0Time>=trgtimes[i]))
3a952328 766 trigInTimeWindow = 1;
767 }
768
769 Double_t tenergy = 0;
770 Double_t tphi=0;
771 Double_t teta=0;
772 for (Int_t i=0;i<3;++i) {
773 Short_t pos = -1;
774 std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
775 if (it!=map.end())
776 pos = it->second;
777 if (pos<0)
b6c599fe 778 continue;
3a952328 779 if (trigInTimeWindow)
5fe1ca23 780 fAmpInTrigger[pos] = amp;
3a952328 781 Float_t eta=-1, phi=-1;
782 fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
783 Double_t en= cells->GetAmplitude(pos);
784 tenergy+=en;
785 teta+=eta*en;
786 tphi+=phi*en;
787 }
788 teta/=tenergy;
789 tphi/=tenergy;
790 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
791 trignew->fE = tenergy;
792 trignew->fEta = teta;
793 trignew->fPhi = tphi;
794 trignew->fAmp = amp;
795 trignew->fMinTime = mintime;
796 trignew->fMaxTime = maxtime;
797 }
798}
799
800//________________________________________________________________________
801void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
802{
803 // Calculate cluster properties
804
805 fClusters->Clear();
806
807 AliVCaloCells *cells = fEsdCells;
808 if (!cells)
809 cells = fAodCells;
810 if (!cells)
811 return;
812
813 TObjArray *clusters = fEsdClusters;
814 if (!clusters)
815 clusters = fAodClusters;
816 if (!clusters)
817 return;
818
819 Int_t ncells = cells->GetNumberOfCells();
820 Int_t nclus = clusters->GetEntries();
821 Int_t ntrks = fSelTracks->GetEntries();
822 Bool_t btracks[6][ntrks];
823 memset(btracks,0,sizeof(btracks));
824
825 std::map<Short_t,Short_t> map;
826 for (Short_t pos=0;pos<ncells;++pos) {
827 Short_t id = cells->GetCellNumber(pos);
828 map[id]=pos;
829 }
830
831 for(Int_t i=0, ncl=0; i<nclus; ++i) {
832 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
833 if (!clus)
834 continue;
835 if (!clus->IsEMCAL())
836 continue;
837 if (clus->E()<fMinE)
838 continue;
839
840 Float_t clsPos[3] = {0};
841 clus->GetPosition(clsPos);
842 TVector3 clsVec(clsPos);
843 Double_t vertex[3] = {0};
844 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
845 TLorentzVector clusterVec;
846 clus->GetMomentum(clusterVec,vertex);
847 Double_t clsEta = clusterVec.Eta();
848
849 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
850 cl->fE = clus->E();
851 cl->fR = clsVec.Perp();
852 cl->fEta = clsVec.Eta();
853 cl->fPhi = clsVec.Phi();
854 cl->fN = clus->GetNCells();
855 cl->fN1 = GetNCells(clus,0.100);
856 cl->fN3 = GetNCells(clus,0.300);
857 Short_t id = -1;
858 Double_t emax = GetMaxCellEnergy(clus, id);
859 cl->fIdMax = id;
860 cl->fEmax = emax;
861 if (clus->GetDistanceToBadChannel()<10000)
862 cl->fDbc = clus->GetDistanceToBadChannel();
863 if (!TMath::IsNaN(clus->GetDispersion()))
864 cl->fDisp = clus->GetDispersion();
865 if (!TMath::IsNaN(clus->GetM20()))
0fbe8d4f 866 cl->fM20 = clus->GetM20();
3a952328 867 if (!TMath::IsNaN(clus->GetM02()))
0fbe8d4f 868 cl->fM02 = clus->GetM02();
3a952328 869 Double_t maxAxis = 0, minAxis = 0;
870 GetSigma(clus,maxAxis,minAxis);
871 clus->SetTOF(maxAxis); // store sigma in TOF
872 cl->fSig = maxAxis;
873 Double_t clusterEcc = 0;
874 if (maxAxis > 0)
875 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
876 clus->SetChi2(clusterEcc); // store ecc in chi2
877 cl->fEcc = clusterEcc;
878 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
879 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
880 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
881 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
882 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
0fbe8d4f 883
3a952328 884 Double_t trigpen = 0;
885 Double_t trignen = 0;
886 for(Int_t j=0; j<cl->fN; ++j) {
887 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
888 Short_t pos = -1;
889 std::map<Short_t,Short_t>::iterator it = map.find(cid);
890 if (it!=map.end())
891 pos = it->second;
892 if (pos<0)
b6c599fe 893 continue;
3a952328 894 if (fAmpInTrigger[pos]>0)
895 trigpen += cells->GetAmplitude(pos);
896 else if (fAmpInTrigger[pos]<0)
897 trignen += cells->GetAmplitude(pos);
b6c599fe 898 }
3a952328 899 if (trigpen>0) {
f3582e89 900 cl->fIsTrigM = 1;
901 cl->fTrigE = trigpen;
3a952328 902 }
903 if (trignen>0) {
f3582e89 904 cl->fIsTrigM = 1;
3a952328 905 cl->fTrigMaskE = trignen;
906 }
0fbe8d4f 907 cl->fIsShared = IsShared(clus);
3a952328 908
909 // track matching
910 Double_t mind2 = 1e10;
911 for(Int_t j = 0; j<ntrks; ++j) {
912 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
b6c599fe 913 if (!track)
914 continue;
3a952328 915
916 if (TMath::Abs(clsEta-track->Eta())>0.5)
b6c599fe 917 continue;
3a952328 918
919 TVector3 vec(clsPos);
920 Int_t index = (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
921 if (btracks[index-4][j]) {
b6c599fe 922 continue;
3a952328 923 }
b6c599fe 924
3a952328 925 Float_t tmpR=-1, tmpZ=-1;
926 if (!fDoTrMatGeom) {
927 AliExternalTrackParam *tParam = 0;
928 if (fEsdEv) {
929 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
930 tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
931 } else
932 tParam = new AliExternalTrackParam(track);
933
934 Double_t bfield[3] = {0};
935 track->GetBxByBz(bfield);
936 Double_t alpha = (index+0.5)*20*TMath::DegToRad();
937 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
938 tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
939 Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
940 if (!ret) {
941 btracks[index-4][j]=1;
942 delete tParam;
943 continue;
b6c599fe 944 }
3a952328 945 Double_t trkPos[3] = {0};
946 tParam->GetXYZ(trkPos); //Get the extrapolated global position
947 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
948 tmpZ = clsPos[2]-trkPos[2];
949 delete tParam;
950 } else {
951 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
952 continue;
953 AliExternalTrackParam tParam(track);
954 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
955 continue;
b6c599fe 956 }
3a952328 957
958 Double_t d2 = tmpR;
959 if (mind2>d2) {
960 mind2=d2;
961 cl->fTrDz = tmpZ;
962 cl->fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
963 cl->fTrEp = clus->E()/track->P();
f3582e89 964 cl->fIsTrackM = 1;
3a952328 965 }
966 }
967
f3582e89 968 if (cl->fIsTrackM) {
3a952328 969 fHMatchDr->Fill(cl->fTrDr);
970 fHMatchDz->Fill(cl->fTrDz);
971 fHMatchEp->Fill(cl->fTrEp);
b6c599fe 972 }
973 }
974}
975
976//________________________________________________________________________
977void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
978{
979 // Calculate track properties.
980
981 fSelPrimTracks->Clear();
982
983 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
984 TClonesArray *tracks = 0;
985 if (fEsdEv) {
986 am->LoadBranch("Tracks");
987 TList *l = fEsdEv->GetList();
988 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
989 } else {
990 am->LoadBranch("tracks");
991 TList *l = fAodEv->GetList();
992 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
993 }
994
296ea9b4 995 if (!tracks)
996 return;
997
0ec74551 998 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 999 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1000 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
0ec74551 1001
296ea9b4 1002 if (fEsdEv) {
b6c599fe 1003 fSelPrimTracks->SetOwner(kTRUE);
0ec74551 1004 am->LoadBranch("PrimaryVertex.");
1005 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1006 am->LoadBranch("SPDVertex.");
1007 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1008 am->LoadBranch("Tracks");
1009 const Int_t Ntracks = tracks->GetEntries();
1010 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1011 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1012 if (!track) {
1013 AliWarning(Form("Could not receive track %d\n", iTracks));
1014 continue;
1015 }
1016 if (fTrCuts && !fTrCuts->IsSelected(track))
296ea9b4 1017 continue;
1018 Double_t eta = track->Eta();
1019 if (eta<-1||eta>1)
1020 continue;
0ec74551 1021 if (track->Phi()<phimin||track->Phi()>phimax)
1022 continue;
2e4d8148 1023
0ec74551 1024 AliESDtrack copyt(*track);
1025 Double_t bfield[3];
1026 copyt.GetBxByBz(bfield);
1027 AliExternalTrackParam tParam;
1028 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1029 if (!relate)
1030 continue;
1031 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
2e4d8148 1032
0ec74551 1033 Double_t p[3] = { 0. };
1034 copyt.GetPxPyPz(p);
1035 Double_t pos[3] = { 0. };
1036 copyt.GetXYZ(pos);
1037 Double_t covTr[21] = { 0. };
1038 copyt.GetCovarianceXYZPxPyPz(covTr);
1039 Double_t pid[10] = { 0. };
1040 copyt.GetESDpid(pid);
1041 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1042 copyt.GetLabel(),
1043 p,
1044 kTRUE,
1045 pos,
1046 kFALSE,
1047 covTr,
1048 (Short_t)copyt.GetSign(),
1049 copyt.GetITSClusterMap(),
1050 pid,
1051 0,/*fPrimaryVertex,*/
1052 kTRUE, // check if this is right
1053 vtx->UsesTrack(copyt.GetID()));
1054 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1055 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1056 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1057 if(ndf>0)
1058 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1059 else
1060 aTrack->SetChi2perNDF(-1);
1061 aTrack->SetFlags(copyt.GetStatus());
1062 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
b6c599fe 1063 fSelPrimTracks->Add(aTrack);
296ea9b4 1064 }
1065 } else {
1066 Int_t ntracks = tracks->GetEntries();
1067 for (Int_t i=0; i<ntracks; ++i) {
1068 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1069 if (!track)
1070 continue;
1071 Double_t eta = track->Eta();
1072 if (eta<-1||eta>1)
1073 continue;
0ec74551 1074 if (track->Phi()<phimin||track->Phi()>phimax)
1075 continue;
b6c599fe 1076 if(track->GetTPCNcls()<fMinNClusPerTr)
296ea9b4 1077 continue;
b6c599fe 1078 //todo: Learn how to set/filter AODs for prim/sec tracks
1079 fSelPrimTracks->Add(track);
296ea9b4 1080 }
1081 }
1082}
1083
1084//________________________________________________________________________
3a952328 1085void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
296ea9b4 1086{
3a952328 1087 // Calculate track properties (including secondaries).
b6c599fe 1088
3a952328 1089 fSelTracks->Clear();
296ea9b4 1090
3a952328 1091 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1092 TClonesArray *tracks = 0;
1093 if (fEsdEv) {
1094 am->LoadBranch("Tracks");
1095 TList *l = fEsdEv->GetList();
1096 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1097 } else {
1098 am->LoadBranch("tracks");
1099 TList *l = fAodEv->GetList();
1100 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1101 }
296ea9b4 1102
3a952328 1103 if (!tracks)
1104 return;
296ea9b4 1105
3a952328 1106 if (fEsdEv) {
1107 const Int_t Ntracks = tracks->GetEntries();
1108 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1109 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1110 if (!track) {
1111 AliWarning(Form("Could not receive track %d\n", iTracks));
1112 continue;
1113 }
1114 if (fTrCuts && !fTrCuts->IsSelected(track))
1115 continue;
1116 Double_t eta = track->Eta();
1117 if (eta<-1||eta>1)
1118 continue;
1119 fSelTracks->Add(track);
1120 }
1121 } else {
1122 Int_t ntracks = tracks->GetEntries();
1123 for (Int_t i=0; i<ntracks; ++i) {
1124 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
296ea9b4 1125 if (!track)
1126 continue;
3a952328 1127 Double_t eta = track->Eta();
1128 if (eta<-1||eta>1)
c4236a58 1129 continue;
3a952328 1130 if(track->GetTPCNcls()<fMinNClusPerTr)
b6c599fe 1131 continue;
2e4d8148 1132
3a952328 1133 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
2e4d8148 1134 AliExternalTrackParam tParam(track);
3a952328 1135 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1136 Double_t trkPos[3];
1137 tParam.GetXYZ(trkPos);
1138 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1139 }
b6c599fe 1140 }
3a952328 1141 fSelTracks->Add(track);
c4236a58 1142 }
296ea9b4 1143 }
1144}
1145
323834f0 1146//________________________________________________________________________
3a952328 1147void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
323834f0 1148{
296ea9b4 1149 // Run custer reconstruction afterburner.
323834f0 1150
1151 AliVCaloCells *cells = fEsdCells;
1152 if (!cells)
1153 cells = fAodCells;
1154
1155 if (!cells)
1156 return;
1157
1158 Int_t ncells = cells->GetNumberOfCells();
1159 if (ncells<=0)
1160 return;
1161
1162 Double_t cellMeanE = 0, cellSigE = 0;
1163 for (Int_t i = 0; i<ncells; ++i) {
1164 Double_t cellE = cells->GetAmplitude(i);
1165 cellMeanE += cellE;
1166 cellSigE += cellE*cellE;
1167 }
1168 cellMeanE /= ncells;
1169 cellSigE /= ncells;
1170 cellSigE -= cellMeanE*cellMeanE;
1171 if (cellSigE<0)
1172 cellSigE = 0;
1173 cellSigE = TMath::Sqrt(cellSigE / ncells);
1174
1175 Double_t subE = cellMeanE - 7*cellSigE;
1176 if (subE<0)
1177 return;
1178
1179 for (Int_t i = 0; i<ncells; ++i) {
1180 Short_t id=-1;
1181 Double_t amp=0,time=0;
1182 if (!cells->GetCell(i, id, amp, time))
1183 continue;
1184 amp -= cellMeanE;
1185 if (amp<0.001)
1186 amp = 0;
1187 cells->SetCell(i, id, amp, time);
1188 }
1189
1190 TObjArray *clusters = fEsdClusters;
1191 if (!clusters)
1192 clusters = fAodClusters;
323834f0 1193 if (!clusters)
1194 return;
1195
1196 Int_t nclus = clusters->GetEntries();
1197 for (Int_t i = 0; i<nclus; ++i) {
1198 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1199 if (!clus->IsEMCAL())
1200 continue;
1201 Int_t nc = clus->GetNCells();
1202 Double_t clusE = 0;
1203 UShort_t ids[100] = {0};
1204 Double_t fra[100] = {0};
1205 for (Int_t j = 0; j<nc; ++j) {
1206 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1207 Double_t cen = cells->GetCellAmplitude(id);
1208 clusE += cen;
1209 if (cen>0) {
1210 ids[nc] = id;
1211 ++nc;
1212 }
1213 }
1214 if (clusE<=0) {
1215 clusters->RemoveAt(i);
1216 continue;
1217 }
1218
1219 for (Int_t j = 0; j<nc; ++j) {
1220 Short_t id = ids[j];
1221 Double_t cen = cells->GetCellAmplitude(id);
1222 fra[j] = cen/clusE;
1223 }
1224 clus->SetE(clusE);
1225 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1226 if (aodclus) {
1227 aodclus->Clear("");
1228 aodclus->SetNCells(nc);
1229 aodclus->SetCellsAmplitudeFraction(fra);
1230 aodclus->SetCellsAbsId(ids);
1231 }
1232 }
1233 clusters->Compress();
1234}
1235
286b47a5 1236//________________________________________________________________________
1237void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1238{
1239 // Fill histograms related to cell properties.
d595acbb 1240
90d5b88b 1241 AliVCaloCells *cells = fEsdCells;
1242 if (!cells)
1243 cells = fAodCells;
1244
296ea9b4 1245 if (!cells)
1246 return;
90d5b88b 1247
296ea9b4 1248 Int_t cellModCount[12] = {0};
1249 Double_t cellMaxE = 0;
1250 Double_t cellMeanE = 0;
1251 Int_t ncells = cells->GetNumberOfCells();
1252 if (ncells==0)
1253 return;
1254
1255 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1256
1257 for (Int_t i = 0; i<ncells; ++i) {
1258 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1259 Double_t cellE = cells->GetAmplitude(i);
1260 fHCellE->Fill(cellE);
1261 if (cellE>cellMaxE)
1262 cellMaxE = cellE;
1263 cellMeanE += cellE;
1264
1265 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1266 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1267 if (!ret) {
1268 AliError(Form("Could not get cell index for %d", absID));
1269 continue;
1270 }
1271 ++cellModCount[iSM];
1272 Int_t iPhi=-1, iEta=-1;
1273 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1274 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1275 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1276 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1277 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1278 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1279 if (fHCellCheckE && fHCellCheckE[absID])
1280 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1281 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1282 }
1283 fHCellH->Fill(cellMaxE);
1284 cellMeanE /= ncells;
1285 fHCellM->Fill(cellMeanE);
1286 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1287 for (Int_t i=0; i<nsm; ++i)
1288 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1289}
1290
1291//________________________________________________________________________
1292void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1293{
90d5b88b 1294 // Fill histograms related to cluster properties.
fa443410 1295
90d5b88b 1296 TObjArray *clusters = fEsdClusters;
1297 if (!clusters)
1298 clusters = fAodClusters;
296ea9b4 1299 if (!clusters)
1300 return;
90d5b88b 1301
296ea9b4 1302 Double_t vertex[3] = {0};
1303 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1304
1305 Int_t nclus = clusters->GetEntries();
1306 for(Int_t i = 0; i<nclus; ++i) {
1307 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1308 if (!clus)
1309 continue;
1310 if (!clus->IsEMCAL())
1311 continue;
90d5b88b 1312 TLorentzVector clusterVec;
296ea9b4 1313 clus->GetMomentum(clusterVec,vertex);
3a952328 1314 Double_t maxAxis = clus->GetTOF(); //sigma
1315 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1316 fHClustEccentricity->Fill(clusterEcc);
1317 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1318 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1319 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1320 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1321 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
788ca675 1322 }
1323}
b3ee6797 1324
788ca675 1325//________________________________________________________________________
1326void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1327{
1328 // Fill ntuple.
1329
1330 if (!fNtuple)
1331 return;
1332
1333 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1334 if (fAodEv) {
1335 fHeader->fRun = fAodEv->GetRunNumber();
1336 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1337 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1338 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1339 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1340 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1341 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1342 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1343 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1344 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1345 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1346 } else {
1347 fHeader->fRun = fEsdEv->GetRunNumber();
1348 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1349 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1350 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1351 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1352 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1353 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1354 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1355 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1356 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1357 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1358 Float_t v0CorrR = 0;
1359 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1360 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1361 if (mult)
1362 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1363 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
788ca675 1364 }
1365 AliCentrality *cent = InputEvent()->GetCentrality();
1366 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1367 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1368 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1369 fHeader->fCqual = cent->GetQuality();
b6c599fe 1370
1371 AliEventplane *ep = InputEvent()->GetEventplane();
1372 if (ep) {
1373 if (ep->GetQVector())
1374 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1375 else
1376 fHeader->fPsi = -1;
b6c599fe 1377 if (ep->GetQsub1()&&ep->GetQsub2())
1378 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1379 else
1380 fHeader->fPsiRes = 0;
1381 }
1382
788ca675 1383 Double_t val = 0;
1384 TString trgclasses(fHeader->fFiredTriggers);
1385 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1386 const char *name = fTrClassNamesArr->At(j)->GetName();
1387 if (trgclasses.Contains(name))
1388 val += TMath::Power(2,j);
1389 }
1390 fHeader->fTcls = (UInt_t)val;
1391
5fe1ca23 1392 fHeader->fNSelTr = fSelTracks->GetEntries();
1393 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
1394
1395 fHeader->fNCells = 0;
1396 fHeader->fNCells1 = 0;
1397 fHeader->fNCells2 = 0;
1398 fHeader->fNCells5 = 0;
1399 fHeader->fMaxCellE = 0;
1400
1401 AliVCaloCells *cells = fEsdCells;
1402 if (!cells)
1403 cells = fAodCells;
1404
1405 if (cells) {
1406 Int_t ncells = cells->GetNumberOfCells();
1407 for(Int_t j=0; j<ncells; ++j) {
1408 Double_t cellen = cells->GetAmplitude(j);
1409 if (cellen>1)
1410 ++fHeader->fNCells1;
1411 if (cellen>2)
1412 ++fHeader->fNCells2;
1413 if (cellen>5)
1414 ++fHeader->fNCells5;
1415 if (cellen>fHeader->fMaxCellE)
1416 fHeader->fMaxCellE = cellen;
1417 }
1418 fHeader->fNCells = ncells;
1419 }
1420
1421 fHeader->fNClus = 0;
1422 fHeader->fNClus1 = 0;
1423 fHeader->fNClus2 = 0;
1424 fHeader->fNClus5 = 0;
1425 fHeader->fMaxClusE = 0;
1426
1427 TObjArray *clusters = fEsdClusters;
1428 if (!clusters)
1429 clusters = fAodClusters;
1430
1431 if (clusters) {
1432 Int_t nclus = clusters->GetEntries();
1433 for(Int_t j=0; j<nclus; ++j) {
1434 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
1435 Double_t clusen = clus->E();
1436 if (clusen>1)
1437 ++fHeader->fNClus1;
1438 if (clusen>2)
1439 ++fHeader->fNClus2;
1440 if (clusen>5)
1441 ++fHeader->fNClus5;
1442 if (clusen>fHeader->fMaxClusE)
1443 fHeader->fMaxClusE = clusen;
1444 }
1445 fHeader->fNClus = nclus;
1446 }
1447
1448
788ca675 1449 if (fAodEv) {
1450 am->LoadBranch("vertices");
1451 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1452 FillVertex(fPrimVert, pv);
1453 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1454 FillVertex(fSpdVert, sv);
1455 } else {
1456 am->LoadBranch("PrimaryVertex.");
1457 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1458 FillVertex(fPrimVert, pv);
1459 am->LoadBranch("SPDVertex.");
1460 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1461 FillVertex(fSpdVert, sv);
1462 am->LoadBranch("TPCVertex.");
1463 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1464 FillVertex(fTpcVert, tv);
fa443410 1465 }
788ca675 1466
788ca675 1467 fNtuple->Fill();
286b47a5 1468}
1469
1470//________________________________________________________________________
1471void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1472{
1473 // Fill histograms related to pions.
286b47a5 1474
296ea9b4 1475 Double_t vertex[3] = {0};
fa443410 1476 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1477
90d5b88b 1478 TObjArray *clusters = fEsdClusters;
1479 if (!clusters)
1480 clusters = fAodClusters;
fa443410 1481
90d5b88b 1482 if (clusters) {
1483 TLorentzVector clusterVec1;
1484 TLorentzVector clusterVec2;
1485 TLorentzVector pionVec;
1486
1487 Int_t nclus = clusters->GetEntries();
fa443410 1488 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1489 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1490 if (!clus1)
1491 continue;
1492 if (!clus1->IsEMCAL())
1493 continue;
3c661da5 1494 if (clus1->E()<fMinE)
6eb6260e 1495 continue;
a49742b5 1496 if (clus1->GetNCells()<fNminCells)
1497 continue;
f224d35b 1498 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1499 continue;
3c661da5 1500 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1501 continue;
fa443410 1502 clus1->GetMomentum(clusterVec1,vertex);
1503 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1504 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1505 if (!clus2)
1506 continue;
1507 if (!clus2->IsEMCAL())
1508 continue;
3c661da5 1509 if (clus2->E()<fMinE)
6eb6260e 1510 continue;
a49742b5 1511 if (clus2->GetNCells()<fNminCells)
1512 continue;
f224d35b 1513 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1514 continue;
3c661da5 1515 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1516 continue;
fa443410 1517 clus2->GetMomentum(clusterVec2,vertex);
1518 pionVec = clusterVec1 + clusterVec2;
1519 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1520 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1521 if (pionZgg < fAsymMax) {
1522 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1523 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1524 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1525 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1526 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1527 fHPionInvMasses[bin]->Fill(pionVec.M());
1528 }
fa443410 1529 }
1530 }
90d5b88b 1531 }
fa443410 1532}
d595acbb 1533
6eb6260e 1534//________________________________________________________________________
323834f0 1535void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1536{
788ca675 1537 // Fill other histograms.
1538}
1539
1540//__________________________________________________________________________________________________
1541void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1542{
1543 // Fill vertex from ESD vertex info.
1544
1545 v->fVx = esdv->GetX();
1546 v->fVy = esdv->GetY();
1547 v->fVz = esdv->GetZ();
1548 v->fVc = esdv->GetNContributors();
1549 v->fDisp = esdv->GetDispersion();
1550 v->fZres = esdv->GetZRes();
1551 v->fChi2 = esdv->GetChi2();
1552 v->fSt = esdv->GetStatus();
1553 v->fIs3D = esdv->IsFromVertexer3D();
1554 v->fIsZ = esdv->IsFromVertexerZ();
1555}
1556
1557//__________________________________________________________________________________________________
1558void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1559{
1560 // Fill vertex from AOD vertex info.
1561
1562 v->fVx = aodv->GetX();
1563 v->fVy = aodv->GetY();
1564 v->fVz = aodv->GetZ();
1565 v->fVc = aodv->GetNContributors();
1566 v->fChi2 = aodv->GetChi2();
6eb6260e 1567}
1568
d595acbb 1569//________________________________________________________________________
296ea9b4 1570Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1571{
1572 // Compute isolation based on cell content.
1573
1574 AliVCaloCells *cells = fEsdCells;
1575 if (!cells)
1576 cells = fAodCells;
1577 if (!cells)
1578 return 0;
1579
1580 Double_t cellIsolation = 0;
1581 Double_t rad2 = radius*radius;
1582 Int_t ncells = cells->GetNumberOfCells();
1583 for (Int_t i = 0; i<ncells; ++i) {
1584 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1585 Float_t eta=-1, phi=-1;
296ea9b4 1586 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1587 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1588 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1589 if(dist>rad2)
1590 continue;
0fbe8d4f 1591 Double_t cellE = cells->GetAmplitude(i);
296ea9b4 1592 cellIsolation += cellE;
1593 }
1594 return cellIsolation;
1595}
1596
1597//________________________________________________________________________
0fbe8d4f 1598Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
1599{
1600 // Get maximum energy of attached cell.
1601
1602 Double_t ret = 0;
1603 Int_t ncells = cluster->GetNCells();
1604 if (fEsdCells) {
1605 for (Int_t i=0; i<ncells; i++) {
1606 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1607 ret += e;
1608 }
1609 } else {
1610 for (Int_t i=0; i<ncells; i++) {
1611 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1612 ret += e;
1613 }
1614 }
1615 return ret;
1616}
1617
1618//________________________________________________________________________
1619Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 1620{
90d5b88b 1621 // Get maximum energy of attached cell.
1622
788ca675 1623 id = -1;
90d5b88b 1624 Double_t maxe = 0;
90d5b88b 1625 Int_t ncells = cluster->GetNCells();
1626 if (fEsdCells) {
1627 for (Int_t i=0; i<ncells; i++) {
27422847 1628 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 1629 if (e>maxe) {
90d5b88b 1630 maxe = e;
788ca675 1631 id = cluster->GetCellAbsId(i);
f0897c18 1632 }
90d5b88b 1633 }
1634 } else {
1635 for (Int_t i=0; i<ncells; i++) {
27422847 1636 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 1637 if (e>maxe)
1638 maxe = e;
788ca675 1639 id = cluster->GetCellAbsId(i);
90d5b88b 1640 }
1641 }
6eb6260e 1642 return maxe;
d595acbb 1643}
1644
1645//________________________________________________________________________
0fbe8d4f 1646void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 1647{
90d5b88b 1648 // Calculate the (E) weighted variance along the longer (eigen) axis.
1649
6bf90832 1650 sigmaMax = 0; // cluster variance along its longer axis
1651 sigmaMin = 0; // cluster variance along its shorter axis
1652 Double_t Ec = c->E(); // cluster energy
296ea9b4 1653 if(Ec<=0)
1654 return;
6bf90832 1655 Double_t Xc = 0 ; // cluster first moment along X
1656 Double_t Yc = 0 ; // cluster first moment along Y
1657 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
1658 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
1659 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 1660
1661 AliVCaloCells *cells = fEsdCells;
1662 if (!cells)
1663 cells = fAodCells;
1664
6eb6260e 1665 if (!cells)
1666 return;
1667
6bf90832 1668 Int_t ncells = c->GetNCells();
6eb6260e 1669 if (ncells==1)
1670 return;
1671
6eb6260e 1672 for(Int_t j=0; j<ncells; ++j) {
6bf90832 1673 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 1674 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 1675 TVector3 pos;
1676 fGeom->GetGlobal(id,pos);
6eb6260e 1677 Xc += cellen*pos.X();
1678 Yc += cellen*pos.Y();
1679 Sxx += cellen*pos.X()*pos.X();
1680 Syy += cellen*pos.Y()*pos.Y();
1681 Sxy += cellen*pos.X()*pos.Y();
1682 }
1683 Xc /= Ec;
1684 Yc /= Ec;
1685 Sxx /= Ec;
1686 Syy /= Ec;
1687 Sxy /= Ec;
1688 Sxx -= Xc*Xc;
1689 Syy -= Yc*Yc;
1690 Sxy -= Xc*Yc;
f0897c18 1691 Sxx = TMath::Abs(Sxx);
1692 Syy = TMath::Abs(Syy);
296ea9b4 1693 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1694 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1695 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1696 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 1697}
90d5b88b 1698
6bf90832 1699//________________________________________________________________________
0fbe8d4f 1700Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 1701{
1702 // Calculate number of attached cells above emin.
1703
6bf90832 1704 AliVCaloCells *cells = fEsdCells;
1705 if (!cells)
1706 cells = fAodCells;
6bf90832 1707 if (!cells)
296ea9b4 1708 return 0;
6bf90832 1709
296ea9b4 1710 Int_t n = 0;
6bf90832 1711 Int_t ncells = c->GetNCells();
1712 for(Int_t j=0; j<ncells; ++j) {
1713 Int_t id = TMath::Abs(c->GetCellAbsId(j));
1714 Double_t cellen = cells->GetCellAmplitude(id);
1715 if (cellen>=emin)
1716 ++n;
1717 }
1718 return n;
1719}
1720
296ea9b4 1721//________________________________________________________________________
b6c599fe 1722Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 1723{
1724 // Compute isolation based on tracks.
1725
1726 Double_t trkIsolation = 0;
1727 Double_t rad2 = radius*radius;
b6c599fe 1728 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 1729 for(Int_t j = 0; j<ntrks; ++j) {
1730 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1731 if (!track)
1732 continue;
b6c599fe 1733 if (track->Pt()<pt)
1734 continue;
296ea9b4 1735 Float_t eta = track->Eta();
1736 Float_t phi = track->Phi();
0fbe8d4f 1737 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1738 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1739 if(dist>rad2)
1740 continue;
1741 trkIsolation += track->Pt();
1742 }
1743 return trkIsolation;
1744}
3a952328 1745
0fbe8d4f 1746//________________________________________________________________________
1747Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
1748{
1749 // Returns if cluster shared across super modules.
1750
1751 AliVCaloCells *cells = fEsdCells;
1752 if (!cells)
1753 cells = fAodCells;
1754 if (!cells)
1755 return 0;
1756
1757 Int_t n = -1;
1758 Int_t ncells = c->GetNCells();
1759 for(Int_t j=0; j<ncells; ++j) {
1760 Int_t id = TMath::Abs(c->GetCellAbsId(j));
1761 Int_t got = id / (24*48);
1762 if (n==-1) {
1763 n = got;
1764 continue;
1765 }
1766 if (got!=n)
1767 return 1;
1768 }
1769 return 0;
1770}
3a952328 1771