]> 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
check if emcal is active
[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.");
44310bce 432 UInt_t mask = fEsdEv->GetESDRun()->GetDetectorsInReco();
433 if ((mask >> 18) & 0x1 == 0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
434 AliError(Form("EMCAL not reconstructed: %u (%u)", mask, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
435 return;
436 }
437 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
43cfaa06 438 } else {
439 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
408dc04c 440 if (!fAodEv) {
441 AliFatal("Neither ESD nor AOD event found");
442 return;
443 }
43cfaa06 444 am->LoadBranch("header");
44310bce 445 offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
446 }
447 if (offtrigger & AliVEvent::kFastOnly) {
448 AliWarning(Form("EMCAL not in fast only partition"));
449 return;
43cfaa06 450 }
b6c599fe 451 if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
27c2e3d9 452 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
453 AliCDBManager *cdb = AliCDBManager::Instance();
454 if (!cdb->IsDefaultStorageSet())
455 cdb->SetDefaultStorage("raw://");
456 Int_t runno = InputEvent()->GetRunNumber();
457 if (runno != cdb->GetRun())
458 cdb->SetRun(runno);
459 AliGeomManager::LoadGeometry();
460 }
461
462 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
463 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
464 if (fEsdEv) {
465 for (Int_t i=0; i<nsm; ++i)
466 fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
467 } else {
468 for (Int_t i=0; i<nsm; ++i)
469 fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
470 }
471 fIsGeoMatsSet = kTRUE;
472 }
473
474 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
475 if (fEsdEv) {
476 const AliESDRun *erun = fEsdEv->GetESDRun();
477 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
478 erun->GetCurrentDip(),
479 AliMagF::kConvLHC,
480 kFALSE,
481 erun->GetBeamEnergy(),
482 erun->GetBeamType());
483 TGeoGlobalMagField::Instance()->SetField(field);
484 } else {
485 Double_t pol = -1; //polarity
486 Double_t be = -1; //beam energy
487 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
488 Int_t runno = fAodEv->GetRunNumber();
489 if (runno>=136851 && runno<138275) {
490 pol = -1;
491 be = 2760;
492 btype = AliMagF::kBeamTypeAA;
493 } else if (runno>=138275 && runno<=139517) {
494 pol = +1;
495 be = 2760;
496 btype = AliMagF::kBeamTypeAA;
497 } else {
498 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
499 }
500 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
501 }
502 }
503
b3ee6797 504 Int_t cut = 1;
505 fHCuts->Fill(cut++);
506
507 TString trgclasses;
508 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
509 if (h) {
510 trgclasses = fEsdEv->GetFiredTriggerClasses();
511 } else {
512 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
513 if (h2)
514 trgclasses = h2->GetFiredTriggerClasses();
515 }
516 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
517 const char *name = fTrClassNamesArr->At(i)->GetName();
518 if (trgclasses.Contains(name))
519 fHTclsBeforeCuts->Fill(1+i);
520 }
521
522 UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
523 if (res==0)
524 return;
525
fa443410 526 fHCuts->Fill(cut++);
286b47a5 527
717fe7de 528 const AliCentrality *centP = InputEvent()->GetCentrality();
529 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
fa443410 530 fHCent->Fill(cent);
717fe7de 531 if (cent<fCentFrom||cent>fCentTo)
286b47a5 532 return;
43cfaa06 533
fa443410 534 fHCuts->Fill(cut++);
286b47a5 535
43cfaa06 536 if (fUseQualFlag) {
537 if (centP->GetQuality()>0)
538 return;
539 }
286b47a5 540
fa443410 541 fHCentQual->Fill(cent);
542 fHCuts->Fill(cut++);
717fe7de 543
43cfaa06 544 if (fEsdEv) {
fa443410 545 am->LoadBranch("PrimaryVertex.");
546 am->LoadBranch("SPDVertex.");
547 am->LoadBranch("TPCVertex.");
43cfaa06 548 } else {
549 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
550 am->LoadBranch("vertices");
3c661da5 551 if (!fAodEv) return;
43cfaa06 552 }
553
554 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
717fe7de 555 if (!vertex)
556 return;
557
fa443410 558 fHVertexZ->Fill(vertex->GetZ());
286b47a5 559
717fe7de 560 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
561 return;
286b47a5 562
fa443410 563 fHCuts->Fill(cut++);
76332037 564 fHVertexZ2->Fill(vertex->GetZ());
717fe7de 565
b3ee6797 566 // count number of accepted events
567 ++fNEvs;
568
569 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
570 const char *name = fTrClassNamesArr->At(i)->GetName();
571 if (trgclasses.Contains(name))
572 fHTclsAfterCuts->Fill(1+i);
573 }
574
43cfaa06 575 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
b6c599fe 576 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 577 fEsdCells = 0; // will be set if ESD input used
b6c599fe 578 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 579 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
580 fAodCells = 0; // will be set if AOD input used
717fe7de 581
43cfaa06 582 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
583 Bool_t clusattached = 0;
584 Bool_t recalibrated = 0;
585 if (1 && !fClusName.IsNull()) {
586 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
587 TObjArray *ts = am->GetTasks();
588 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
589 if (cltask && cltask->GetClusters()) {
d595acbb 590 fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
43cfaa06 591 clusattached = cltask->GetAttachClusters();
592 if (cltask->GetCalibData()!=0)
593 recalibrated = kTRUE;
594 }
595 }
596 if (1 && AODEvent() && !fClusName.IsNull()) {
717fe7de 597 TList *l = AODEvent()->GetList();
598 TClonesArray *clus = 0;
599 if (l) {
600 clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
601 fAodClusters = clus;
602 }
ea3fd2d5 603 }
717fe7de 604
605 if (fEsdEv) { // ESD input mode
43cfaa06 606 if (1 && (!fRecPoints||clusattached)) {
607 if (!clusattached)
608 am->LoadBranch("CaloClusters");
717fe7de 609 TList *l = fEsdEv->GetList();
610 TClonesArray *clus = 0;
611 if (l) {
612 clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
613 fEsdClusters = clus;
ea3fd2d5 614 }
615 }
43cfaa06 616 if (1) {
617 if (!recalibrated)
618 am->LoadBranch("EMCALCells.");
717fe7de 619 fEsdCells = fEsdEv->GetEMCALCells();
620 }
621 } else if (fAodEv) { // AOD input mode
43cfaa06 622 if (1 && (!fAodClusters || clusattached)) {
623 if (!clusattached)
624 am->LoadBranch("caloClusters");
717fe7de 625 TList *l = fAodEv->GetList();
626 TClonesArray *clus = 0;
627 if (l) {
628 clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
629 fAodClusters = clus;
ea3fd2d5 630 }
631 }
43cfaa06 632 if (1) {
633 if (!recalibrated)
634 am->LoadBranch("emcalCells");
717fe7de 635 fAodCells = fAodEv->GetEMCALCells();
636 }
637 } else {
638 AliFatal("Impossible to not have either pointer to ESD or AOD event");
639 }
ea3fd2d5 640
43cfaa06 641 if (1) {
642 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
643 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
644 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
645 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
646 AliDebug(2,Form("fAodCells set: %p", fAodCells));
647 }
648
a49742b5 649 if (fDoAfterburner)
650 ClusterAfterburner();
6eb6260e 651
3a952328 652 CalcCaloTriggers();
b6c599fe 653 CalcPrimTracks();
3a952328 654 CalcTracks();
296ea9b4 655 CalcClusterProps();
656
286b47a5 657 FillCellHists();
3a952328 658 if (!fTrainMode) {
659 FillClusHists();
660 FillPionHists();
661 FillOtherHists();
662 }
788ca675 663 FillNtuple();
ea3fd2d5 664
3a952328 665 if (fTrainMode) {
666 fSelTracks->Clear();
667 fSelPrimTracks->Clear();
668 fTriggers->Clear();
669 }
de4cee41 670
717fe7de 671 PostData(1, fOutput);
ea3fd2d5 672}
673
674//________________________________________________________________________
675void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
676{
717fe7de 677 // Terminate called at the end of analysis.
f5d4ab70 678
679 if (fNtuple) {
680 TFile *f = OpenFile(1);
681 if (f)
682 fNtuple->Write();
683 }
d9f26424 684
b3ee6797 685 AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
ea3fd2d5 686}
ea3fd2d5 687
296ea9b4 688//________________________________________________________________________
3a952328 689void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
296ea9b4 690{
3a952328 691 // Calculate triggers
296ea9b4 692
3a952328 693 fTriggers->Clear();
296ea9b4 694
3a952328 695 AliVCaloCells *cells = fEsdCells;
696 if (!cells)
697 cells = fAodCells;
698 if (!cells)
699 return;
700
701 Int_t ncells = cells->GetNumberOfCells();
702 if (ncells<=0)
703 return;
704
705 if (ncells>fNAmpInTrigger) {
706 delete [] fAmpInTrigger;
707 fAmpInTrigger = new Float_t[ncells];
708 fNAmpInTrigger = ncells;
296ea9b4 709 }
3a952328 710 for (Int_t i=0;i<ncells;++i)
711 fAmpInTrigger[i] = 0;
296ea9b4 712
3a952328 713 std::map<Short_t,Short_t> map;
714 for (Short_t pos=0;pos<ncells;++pos) {
715 Short_t id = cells->GetCellNumber(pos);
716 map[id]=pos;
717 }
718
719 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
720 am->LoadBranch("EMCALTrigger.");
721
722 AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
723 if (!triggers)
724 return;
725 if (triggers->GetEntries()<=0)
b6c599fe 726 return;
727
3a952328 728 triggers->Reset();
729 Int_t ntrigs=0;
730 while (triggers->Next()) {
731 Int_t gCol=0, gRow=0, ntimes=0;
732 triggers->GetPosition(gCol,gRow);
733 triggers->GetNL0Times(ntimes);
734 if (ntimes<1)
735 continue;
736 Float_t amp=0;
737 triggers->GetAmplitude(amp);
738 Int_t find = -1;
739 fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
740 if (find<0)
741 continue;
742 Int_t cidx[4] = {-1};
743 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
744 if (!ret)
745 continue;
746 Int_t trgtimes[25];
747 triggers->GetL0Times(trgtimes);
748 Int_t mintime = trgtimes[0];
749 Int_t maxtime = trgtimes[0];
750 Bool_t trigInTimeWindow = 0;
751 for (Int_t i=0;i<ntimes;++i) {
752 if (trgtimes[i]<mintime)
753 mintime = trgtimes[i];
754 if (maxtime<trgtimes[i])
755 maxtime = trgtimes[i];
f2b8fca6 756 if ((fMinL0Time<=trgtimes[i]) && (fMaxL0Time>=trgtimes[i]))
3a952328 757 trigInTimeWindow = 1;
758 }
759
760 Double_t tenergy = 0;
761 Double_t tphi=0;
762 Double_t teta=0;
763 for (Int_t i=0;i<3;++i) {
764 Short_t pos = -1;
765 std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
766 if (it!=map.end())
767 pos = it->second;
768 if (pos<0)
b6c599fe 769 continue;
3a952328 770 if (trigInTimeWindow)
5fe1ca23 771 fAmpInTrigger[pos] = amp;
3a952328 772 Float_t eta=-1, phi=-1;
773 fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
774 Double_t en= cells->GetAmplitude(pos);
775 tenergy+=en;
776 teta+=eta*en;
777 tphi+=phi*en;
778 }
779 teta/=tenergy;
780 tphi/=tenergy;
781 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
782 trignew->fE = tenergy;
783 trignew->fEta = teta;
784 trignew->fPhi = tphi;
785 trignew->fAmp = amp;
786 trignew->fMinTime = mintime;
787 trignew->fMaxTime = maxtime;
788 }
789}
790
791//________________________________________________________________________
792void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
793{
794 // Calculate cluster properties
795
796 fClusters->Clear();
797
798 AliVCaloCells *cells = fEsdCells;
799 if (!cells)
800 cells = fAodCells;
801 if (!cells)
802 return;
803
804 TObjArray *clusters = fEsdClusters;
805 if (!clusters)
806 clusters = fAodClusters;
807 if (!clusters)
808 return;
809
810 Int_t ncells = cells->GetNumberOfCells();
811 Int_t nclus = clusters->GetEntries();
812 Int_t ntrks = fSelTracks->GetEntries();
813 Bool_t btracks[6][ntrks];
814 memset(btracks,0,sizeof(btracks));
815
816 std::map<Short_t,Short_t> map;
817 for (Short_t pos=0;pos<ncells;++pos) {
818 Short_t id = cells->GetCellNumber(pos);
819 map[id]=pos;
820 }
821
822 for(Int_t i=0, ncl=0; i<nclus; ++i) {
823 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
824 if (!clus)
825 continue;
826 if (!clus->IsEMCAL())
827 continue;
828 if (clus->E()<fMinE)
829 continue;
830
831 Float_t clsPos[3] = {0};
832 clus->GetPosition(clsPos);
833 TVector3 clsVec(clsPos);
834 Double_t vertex[3] = {0};
835 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
836 TLorentzVector clusterVec;
837 clus->GetMomentum(clusterVec,vertex);
838 Double_t clsEta = clusterVec.Eta();
839
840 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
841 cl->fE = clus->E();
842 cl->fR = clsVec.Perp();
843 cl->fEta = clsVec.Eta();
844 cl->fPhi = clsVec.Phi();
845 cl->fN = clus->GetNCells();
846 cl->fN1 = GetNCells(clus,0.100);
847 cl->fN3 = GetNCells(clus,0.300);
848 Short_t id = -1;
849 Double_t emax = GetMaxCellEnergy(clus, id);
850 cl->fIdMax = id;
851 cl->fEmax = emax;
852 if (clus->GetDistanceToBadChannel()<10000)
853 cl->fDbc = clus->GetDistanceToBadChannel();
854 if (!TMath::IsNaN(clus->GetDispersion()))
855 cl->fDisp = clus->GetDispersion();
856 if (!TMath::IsNaN(clus->GetM20()))
0fbe8d4f 857 cl->fM20 = clus->GetM20();
3a952328 858 if (!TMath::IsNaN(clus->GetM02()))
0fbe8d4f 859 cl->fM02 = clus->GetM02();
3a952328 860 Double_t maxAxis = 0, minAxis = 0;
861 GetSigma(clus,maxAxis,minAxis);
862 clus->SetTOF(maxAxis); // store sigma in TOF
863 cl->fSig = maxAxis;
864 Double_t clusterEcc = 0;
865 if (maxAxis > 0)
866 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
867 clus->SetChi2(clusterEcc); // store ecc in chi2
868 cl->fEcc = clusterEcc;
869 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
870 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
871 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
872 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
873 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
0fbe8d4f 874
3a952328 875 Double_t trigpen = 0;
876 Double_t trignen = 0;
877 for(Int_t j=0; j<cl->fN; ++j) {
878 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
879 Short_t pos = -1;
880 std::map<Short_t,Short_t>::iterator it = map.find(cid);
881 if (it!=map.end())
882 pos = it->second;
883 if (pos<0)
b6c599fe 884 continue;
3a952328 885 if (fAmpInTrigger[pos]>0)
886 trigpen += cells->GetAmplitude(pos);
887 else if (fAmpInTrigger[pos]<0)
888 trignen += cells->GetAmplitude(pos);
b6c599fe 889 }
3a952328 890 if (trigpen>0) {
f3582e89 891 cl->fIsTrigM = 1;
892 cl->fTrigE = trigpen;
3a952328 893 }
894 if (trignen>0) {
f3582e89 895 cl->fIsTrigM = 1;
3a952328 896 cl->fTrigMaskE = trignen;
897 }
0fbe8d4f 898 cl->fIsShared = IsShared(clus);
3a952328 899
900 // track matching
901 Double_t mind2 = 1e10;
902 for(Int_t j = 0; j<ntrks; ++j) {
903 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
b6c599fe 904 if (!track)
905 continue;
3a952328 906
907 if (TMath::Abs(clsEta-track->Eta())>0.5)
b6c599fe 908 continue;
3a952328 909
910 TVector3 vec(clsPos);
911 Int_t index = (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
912 if (btracks[index-4][j]) {
b6c599fe 913 continue;
3a952328 914 }
b6c599fe 915
3a952328 916 Float_t tmpR=-1, tmpZ=-1;
917 if (!fDoTrMatGeom) {
918 AliExternalTrackParam *tParam = 0;
919 if (fEsdEv) {
920 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
921 tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
922 } else
923 tParam = new AliExternalTrackParam(track);
924
925 Double_t bfield[3] = {0};
926 track->GetBxByBz(bfield);
927 Double_t alpha = (index+0.5)*20*TMath::DegToRad();
928 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
929 tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
930 Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
931 if (!ret) {
932 btracks[index-4][j]=1;
933 delete tParam;
934 continue;
b6c599fe 935 }
3a952328 936 Double_t trkPos[3] = {0};
937 tParam->GetXYZ(trkPos); //Get the extrapolated global position
938 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
939 tmpZ = clsPos[2]-trkPos[2];
940 delete tParam;
941 } else {
942 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
943 continue;
944 AliExternalTrackParam tParam(track);
945 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
946 continue;
b6c599fe 947 }
3a952328 948
949 Double_t d2 = tmpR;
950 if (mind2>d2) {
951 mind2=d2;
952 cl->fTrDz = tmpZ;
953 cl->fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
954 cl->fTrEp = clus->E()/track->P();
f3582e89 955 cl->fIsTrackM = 1;
3a952328 956 }
957 }
958
f3582e89 959 if (cl->fIsTrackM) {
3a952328 960 fHMatchDr->Fill(cl->fTrDr);
961 fHMatchDz->Fill(cl->fTrDz);
962 fHMatchEp->Fill(cl->fTrEp);
b6c599fe 963 }
964 }
965}
966
967//________________________________________________________________________
968void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
969{
970 // Calculate track properties.
971
972 fSelPrimTracks->Clear();
973
974 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
975 TClonesArray *tracks = 0;
976 if (fEsdEv) {
977 am->LoadBranch("Tracks");
978 TList *l = fEsdEv->GetList();
979 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
980 } else {
981 am->LoadBranch("tracks");
982 TList *l = fAodEv->GetList();
983 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
984 }
985
296ea9b4 986 if (!tracks)
987 return;
988
0ec74551 989 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 990 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
991 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
0ec74551 992
296ea9b4 993 if (fEsdEv) {
b6c599fe 994 fSelPrimTracks->SetOwner(kTRUE);
0ec74551 995 am->LoadBranch("PrimaryVertex.");
996 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
997 am->LoadBranch("SPDVertex.");
998 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
999 am->LoadBranch("Tracks");
1000 const Int_t Ntracks = tracks->GetEntries();
1001 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1002 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1003 if (!track) {
1004 AliWarning(Form("Could not receive track %d\n", iTracks));
1005 continue;
1006 }
1007 if (fTrCuts && !fTrCuts->IsSelected(track))
296ea9b4 1008 continue;
1009 Double_t eta = track->Eta();
1010 if (eta<-1||eta>1)
1011 continue;
0ec74551 1012 if (track->Phi()<phimin||track->Phi()>phimax)
1013 continue;
2e4d8148 1014
0ec74551 1015 AliESDtrack copyt(*track);
1016 Double_t bfield[3];
1017 copyt.GetBxByBz(bfield);
1018 AliExternalTrackParam tParam;
1019 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1020 if (!relate)
1021 continue;
1022 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
2e4d8148 1023
0ec74551 1024 Double_t p[3] = { 0. };
1025 copyt.GetPxPyPz(p);
1026 Double_t pos[3] = { 0. };
1027 copyt.GetXYZ(pos);
1028 Double_t covTr[21] = { 0. };
1029 copyt.GetCovarianceXYZPxPyPz(covTr);
1030 Double_t pid[10] = { 0. };
1031 copyt.GetESDpid(pid);
1032 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1033 copyt.GetLabel(),
1034 p,
1035 kTRUE,
1036 pos,
1037 kFALSE,
1038 covTr,
1039 (Short_t)copyt.GetSign(),
1040 copyt.GetITSClusterMap(),
1041 pid,
1042 0,/*fPrimaryVertex,*/
1043 kTRUE, // check if this is right
1044 vtx->UsesTrack(copyt.GetID()));
1045 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1046 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1047 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1048 if(ndf>0)
1049 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1050 else
1051 aTrack->SetChi2perNDF(-1);
1052 aTrack->SetFlags(copyt.GetStatus());
1053 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
b6c599fe 1054 fSelPrimTracks->Add(aTrack);
296ea9b4 1055 }
1056 } else {
1057 Int_t ntracks = tracks->GetEntries();
1058 for (Int_t i=0; i<ntracks; ++i) {
1059 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1060 if (!track)
1061 continue;
1062 Double_t eta = track->Eta();
1063 if (eta<-1||eta>1)
1064 continue;
0ec74551 1065 if (track->Phi()<phimin||track->Phi()>phimax)
1066 continue;
b6c599fe 1067 if(track->GetTPCNcls()<fMinNClusPerTr)
296ea9b4 1068 continue;
b6c599fe 1069 //todo: Learn how to set/filter AODs for prim/sec tracks
1070 fSelPrimTracks->Add(track);
296ea9b4 1071 }
1072 }
1073}
1074
1075//________________________________________________________________________
3a952328 1076void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
296ea9b4 1077{
3a952328 1078 // Calculate track properties (including secondaries).
b6c599fe 1079
3a952328 1080 fSelTracks->Clear();
296ea9b4 1081
3a952328 1082 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1083 TClonesArray *tracks = 0;
1084 if (fEsdEv) {
1085 am->LoadBranch("Tracks");
1086 TList *l = fEsdEv->GetList();
1087 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1088 } else {
1089 am->LoadBranch("tracks");
1090 TList *l = fAodEv->GetList();
1091 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1092 }
296ea9b4 1093
3a952328 1094 if (!tracks)
1095 return;
296ea9b4 1096
3a952328 1097 if (fEsdEv) {
1098 const Int_t Ntracks = tracks->GetEntries();
1099 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1100 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1101 if (!track) {
1102 AliWarning(Form("Could not receive track %d\n", iTracks));
1103 continue;
1104 }
1105 if (fTrCuts && !fTrCuts->IsSelected(track))
1106 continue;
1107 Double_t eta = track->Eta();
1108 if (eta<-1||eta>1)
1109 continue;
1110 fSelTracks->Add(track);
1111 }
1112 } else {
1113 Int_t ntracks = tracks->GetEntries();
1114 for (Int_t i=0; i<ntracks; ++i) {
1115 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
296ea9b4 1116 if (!track)
1117 continue;
3a952328 1118 Double_t eta = track->Eta();
1119 if (eta<-1||eta>1)
c4236a58 1120 continue;
3a952328 1121 if(track->GetTPCNcls()<fMinNClusPerTr)
b6c599fe 1122 continue;
2e4d8148 1123
3a952328 1124 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
2e4d8148 1125 AliExternalTrackParam tParam(track);
3a952328 1126 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1127 Double_t trkPos[3];
1128 tParam.GetXYZ(trkPos);
1129 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1130 }
b6c599fe 1131 }
3a952328 1132 fSelTracks->Add(track);
c4236a58 1133 }
296ea9b4 1134 }
1135}
1136
323834f0 1137//________________________________________________________________________
3a952328 1138void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
323834f0 1139{
296ea9b4 1140 // Run custer reconstruction afterburner.
323834f0 1141
1142 AliVCaloCells *cells = fEsdCells;
1143 if (!cells)
1144 cells = fAodCells;
1145
1146 if (!cells)
1147 return;
1148
1149 Int_t ncells = cells->GetNumberOfCells();
1150 if (ncells<=0)
1151 return;
1152
1153 Double_t cellMeanE = 0, cellSigE = 0;
1154 for (Int_t i = 0; i<ncells; ++i) {
1155 Double_t cellE = cells->GetAmplitude(i);
1156 cellMeanE += cellE;
1157 cellSigE += cellE*cellE;
1158 }
1159 cellMeanE /= ncells;
1160 cellSigE /= ncells;
1161 cellSigE -= cellMeanE*cellMeanE;
1162 if (cellSigE<0)
1163 cellSigE = 0;
1164 cellSigE = TMath::Sqrt(cellSigE / ncells);
1165
1166 Double_t subE = cellMeanE - 7*cellSigE;
1167 if (subE<0)
1168 return;
1169
1170 for (Int_t i = 0; i<ncells; ++i) {
1171 Short_t id=-1;
1172 Double_t amp=0,time=0;
1173 if (!cells->GetCell(i, id, amp, time))
1174 continue;
1175 amp -= cellMeanE;
1176 if (amp<0.001)
1177 amp = 0;
1178 cells->SetCell(i, id, amp, time);
1179 }
1180
1181 TObjArray *clusters = fEsdClusters;
1182 if (!clusters)
1183 clusters = fAodClusters;
323834f0 1184 if (!clusters)
1185 return;
1186
1187 Int_t nclus = clusters->GetEntries();
1188 for (Int_t i = 0; i<nclus; ++i) {
1189 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1190 if (!clus->IsEMCAL())
1191 continue;
1192 Int_t nc = clus->GetNCells();
1193 Double_t clusE = 0;
1194 UShort_t ids[100] = {0};
1195 Double_t fra[100] = {0};
1196 for (Int_t j = 0; j<nc; ++j) {
1197 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1198 Double_t cen = cells->GetCellAmplitude(id);
1199 clusE += cen;
1200 if (cen>0) {
1201 ids[nc] = id;
1202 ++nc;
1203 }
1204 }
1205 if (clusE<=0) {
1206 clusters->RemoveAt(i);
1207 continue;
1208 }
1209
1210 for (Int_t j = 0; j<nc; ++j) {
1211 Short_t id = ids[j];
1212 Double_t cen = cells->GetCellAmplitude(id);
1213 fra[j] = cen/clusE;
1214 }
1215 clus->SetE(clusE);
1216 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1217 if (aodclus) {
1218 aodclus->Clear("");
1219 aodclus->SetNCells(nc);
1220 aodclus->SetCellsAmplitudeFraction(fra);
1221 aodclus->SetCellsAbsId(ids);
1222 }
1223 }
1224 clusters->Compress();
1225}
1226
286b47a5 1227//________________________________________________________________________
1228void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1229{
1230 // Fill histograms related to cell properties.
d595acbb 1231
90d5b88b 1232 AliVCaloCells *cells = fEsdCells;
1233 if (!cells)
1234 cells = fAodCells;
1235
296ea9b4 1236 if (!cells)
1237 return;
90d5b88b 1238
296ea9b4 1239 Int_t cellModCount[12] = {0};
1240 Double_t cellMaxE = 0;
1241 Double_t cellMeanE = 0;
1242 Int_t ncells = cells->GetNumberOfCells();
1243 if (ncells==0)
1244 return;
1245
1246 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1247
1248 for (Int_t i = 0; i<ncells; ++i) {
1249 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1250 Double_t cellE = cells->GetAmplitude(i);
1251 fHCellE->Fill(cellE);
1252 if (cellE>cellMaxE)
1253 cellMaxE = cellE;
1254 cellMeanE += cellE;
1255
1256 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1257 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1258 if (!ret) {
1259 AliError(Form("Could not get cell index for %d", absID));
1260 continue;
1261 }
1262 ++cellModCount[iSM];
1263 Int_t iPhi=-1, iEta=-1;
1264 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1265 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1266 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1267 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1268 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1269 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1270 if (fHCellCheckE && fHCellCheckE[absID])
1271 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1272 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1273 }
1274 fHCellH->Fill(cellMaxE);
1275 cellMeanE /= ncells;
1276 fHCellM->Fill(cellMeanE);
1277 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1278 for (Int_t i=0; i<nsm; ++i)
1279 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1280}
1281
1282//________________________________________________________________________
1283void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1284{
90d5b88b 1285 // Fill histograms related to cluster properties.
fa443410 1286
90d5b88b 1287 TObjArray *clusters = fEsdClusters;
1288 if (!clusters)
1289 clusters = fAodClusters;
296ea9b4 1290 if (!clusters)
1291 return;
90d5b88b 1292
296ea9b4 1293 Double_t vertex[3] = {0};
1294 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1295
1296 Int_t nclus = clusters->GetEntries();
1297 for(Int_t i = 0; i<nclus; ++i) {
1298 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1299 if (!clus)
1300 continue;
1301 if (!clus->IsEMCAL())
1302 continue;
90d5b88b 1303 TLorentzVector clusterVec;
296ea9b4 1304 clus->GetMomentum(clusterVec,vertex);
3a952328 1305 Double_t maxAxis = clus->GetTOF(); //sigma
1306 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1307 fHClustEccentricity->Fill(clusterEcc);
1308 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1309 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1310 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1311 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1312 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
788ca675 1313 }
1314}
b3ee6797 1315
788ca675 1316//________________________________________________________________________
1317void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1318{
1319 // Fill ntuple.
1320
1321 if (!fNtuple)
1322 return;
1323
1324 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1325 if (fAodEv) {
1326 fHeader->fRun = fAodEv->GetRunNumber();
1327 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1328 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1329 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1330 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1331 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1332 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1333 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1334 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1335 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1336 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1337 } else {
1338 fHeader->fRun = fEsdEv->GetRunNumber();
1339 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1340 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1341 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1342 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1343 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1344 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1345 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1346 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1347 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1348 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1349 Float_t v0CorrR = 0;
1350 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1351 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1352 if (mult)
1353 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1354 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
788ca675 1355 }
1356 AliCentrality *cent = InputEvent()->GetCentrality();
1357 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1358 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1359 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1360 fHeader->fCqual = cent->GetQuality();
b6c599fe 1361
1362 AliEventplane *ep = InputEvent()->GetEventplane();
1363 if (ep) {
1364 if (ep->GetQVector())
1365 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1366 else
1367 fHeader->fPsi = -1;
b6c599fe 1368 if (ep->GetQsub1()&&ep->GetQsub2())
1369 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1370 else
1371 fHeader->fPsiRes = 0;
1372 }
1373
788ca675 1374 Double_t val = 0;
1375 TString trgclasses(fHeader->fFiredTriggers);
1376 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1377 const char *name = fTrClassNamesArr->At(j)->GetName();
1378 if (trgclasses.Contains(name))
1379 val += TMath::Power(2,j);
1380 }
1381 fHeader->fTcls = (UInt_t)val;
1382
5fe1ca23 1383 fHeader->fNSelTr = fSelTracks->GetEntries();
1384 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
1385
1386 fHeader->fNCells = 0;
1387 fHeader->fNCells1 = 0;
1388 fHeader->fNCells2 = 0;
1389 fHeader->fNCells5 = 0;
1390 fHeader->fMaxCellE = 0;
1391
1392 AliVCaloCells *cells = fEsdCells;
1393 if (!cells)
1394 cells = fAodCells;
1395
1396 if (cells) {
1397 Int_t ncells = cells->GetNumberOfCells();
1398 for(Int_t j=0; j<ncells; ++j) {
1399 Double_t cellen = cells->GetAmplitude(j);
1400 if (cellen>1)
1401 ++fHeader->fNCells1;
1402 if (cellen>2)
1403 ++fHeader->fNCells2;
1404 if (cellen>5)
1405 ++fHeader->fNCells5;
1406 if (cellen>fHeader->fMaxCellE)
1407 fHeader->fMaxCellE = cellen;
1408 }
1409 fHeader->fNCells = ncells;
1410 }
1411
1412 fHeader->fNClus = 0;
1413 fHeader->fNClus1 = 0;
1414 fHeader->fNClus2 = 0;
1415 fHeader->fNClus5 = 0;
1416 fHeader->fMaxClusE = 0;
1417
1418 TObjArray *clusters = fEsdClusters;
1419 if (!clusters)
1420 clusters = fAodClusters;
1421
1422 if (clusters) {
1423 Int_t nclus = clusters->GetEntries();
1424 for(Int_t j=0; j<nclus; ++j) {
1425 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
1426 Double_t clusen = clus->E();
1427 if (clusen>1)
1428 ++fHeader->fNClus1;
1429 if (clusen>2)
1430 ++fHeader->fNClus2;
1431 if (clusen>5)
1432 ++fHeader->fNClus5;
1433 if (clusen>fHeader->fMaxClusE)
1434 fHeader->fMaxClusE = clusen;
1435 }
1436 fHeader->fNClus = nclus;
1437 }
1438
1439
788ca675 1440 if (fAodEv) {
1441 am->LoadBranch("vertices");
1442 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1443 FillVertex(fPrimVert, pv);
1444 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1445 FillVertex(fSpdVert, sv);
1446 } else {
1447 am->LoadBranch("PrimaryVertex.");
1448 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1449 FillVertex(fPrimVert, pv);
1450 am->LoadBranch("SPDVertex.");
1451 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1452 FillVertex(fSpdVert, sv);
1453 am->LoadBranch("TPCVertex.");
1454 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1455 FillVertex(fTpcVert, tv);
fa443410 1456 }
788ca675 1457
788ca675 1458 fNtuple->Fill();
286b47a5 1459}
1460
1461//________________________________________________________________________
1462void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1463{
1464 // Fill histograms related to pions.
286b47a5 1465
296ea9b4 1466 Double_t vertex[3] = {0};
fa443410 1467 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1468
90d5b88b 1469 TObjArray *clusters = fEsdClusters;
1470 if (!clusters)
1471 clusters = fAodClusters;
fa443410 1472
90d5b88b 1473 if (clusters) {
1474 TLorentzVector clusterVec1;
1475 TLorentzVector clusterVec2;
1476 TLorentzVector pionVec;
1477
1478 Int_t nclus = clusters->GetEntries();
fa443410 1479 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1480 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1481 if (!clus1)
1482 continue;
1483 if (!clus1->IsEMCAL())
1484 continue;
3c661da5 1485 if (clus1->E()<fMinE)
6eb6260e 1486 continue;
a49742b5 1487 if (clus1->GetNCells()<fNminCells)
1488 continue;
f224d35b 1489 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1490 continue;
3c661da5 1491 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1492 continue;
fa443410 1493 clus1->GetMomentum(clusterVec1,vertex);
1494 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1495 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1496 if (!clus2)
1497 continue;
1498 if (!clus2->IsEMCAL())
1499 continue;
3c661da5 1500 if (clus2->E()<fMinE)
6eb6260e 1501 continue;
a49742b5 1502 if (clus2->GetNCells()<fNminCells)
1503 continue;
f224d35b 1504 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1505 continue;
3c661da5 1506 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1507 continue;
fa443410 1508 clus2->GetMomentum(clusterVec2,vertex);
1509 pionVec = clusterVec1 + clusterVec2;
1510 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1511 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1512 if (pionZgg < fAsymMax) {
1513 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1514 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1515 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1516 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1517 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1518 fHPionInvMasses[bin]->Fill(pionVec.M());
1519 }
fa443410 1520 }
1521 }
90d5b88b 1522 }
fa443410 1523}
d595acbb 1524
6eb6260e 1525//________________________________________________________________________
323834f0 1526void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1527{
788ca675 1528 // Fill other histograms.
1529}
1530
1531//__________________________________________________________________________________________________
1532void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1533{
1534 // Fill vertex from ESD vertex info.
1535
1536 v->fVx = esdv->GetX();
1537 v->fVy = esdv->GetY();
1538 v->fVz = esdv->GetZ();
1539 v->fVc = esdv->GetNContributors();
1540 v->fDisp = esdv->GetDispersion();
1541 v->fZres = esdv->GetZRes();
1542 v->fChi2 = esdv->GetChi2();
1543 v->fSt = esdv->GetStatus();
1544 v->fIs3D = esdv->IsFromVertexer3D();
1545 v->fIsZ = esdv->IsFromVertexerZ();
1546}
1547
1548//__________________________________________________________________________________________________
1549void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1550{
1551 // Fill vertex from AOD vertex info.
1552
1553 v->fVx = aodv->GetX();
1554 v->fVy = aodv->GetY();
1555 v->fVz = aodv->GetZ();
1556 v->fVc = aodv->GetNContributors();
1557 v->fChi2 = aodv->GetChi2();
6eb6260e 1558}
1559
d595acbb 1560//________________________________________________________________________
296ea9b4 1561Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1562{
1563 // Compute isolation based on cell content.
1564
1565 AliVCaloCells *cells = fEsdCells;
1566 if (!cells)
1567 cells = fAodCells;
1568 if (!cells)
1569 return 0;
1570
1571 Double_t cellIsolation = 0;
1572 Double_t rad2 = radius*radius;
1573 Int_t ncells = cells->GetNumberOfCells();
1574 for (Int_t i = 0; i<ncells; ++i) {
1575 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1576 Float_t eta=-1, phi=-1;
296ea9b4 1577 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1578 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1579 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1580 if(dist>rad2)
1581 continue;
0fbe8d4f 1582 Double_t cellE = cells->GetAmplitude(i);
296ea9b4 1583 cellIsolation += cellE;
1584 }
1585 return cellIsolation;
1586}
1587
1588//________________________________________________________________________
0fbe8d4f 1589Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
1590{
1591 // Get maximum energy of attached cell.
1592
1593 Double_t ret = 0;
1594 Int_t ncells = cluster->GetNCells();
1595 if (fEsdCells) {
1596 for (Int_t i=0; i<ncells; i++) {
1597 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1598 ret += e;
1599 }
1600 } else {
1601 for (Int_t i=0; i<ncells; i++) {
1602 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1603 ret += e;
1604 }
1605 }
1606 return ret;
1607}
1608
1609//________________________________________________________________________
1610Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 1611{
90d5b88b 1612 // Get maximum energy of attached cell.
1613
788ca675 1614 id = -1;
90d5b88b 1615 Double_t maxe = 0;
90d5b88b 1616 Int_t ncells = cluster->GetNCells();
1617 if (fEsdCells) {
1618 for (Int_t i=0; i<ncells; i++) {
27422847 1619 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 1620 if (e>maxe) {
90d5b88b 1621 maxe = e;
788ca675 1622 id = cluster->GetCellAbsId(i);
f0897c18 1623 }
90d5b88b 1624 }
1625 } else {
1626 for (Int_t i=0; i<ncells; i++) {
27422847 1627 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 1628 if (e>maxe)
1629 maxe = e;
788ca675 1630 id = cluster->GetCellAbsId(i);
90d5b88b 1631 }
1632 }
6eb6260e 1633 return maxe;
d595acbb 1634}
1635
1636//________________________________________________________________________
0fbe8d4f 1637void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 1638{
90d5b88b 1639 // Calculate the (E) weighted variance along the longer (eigen) axis.
1640
6bf90832 1641 sigmaMax = 0; // cluster variance along its longer axis
1642 sigmaMin = 0; // cluster variance along its shorter axis
1643 Double_t Ec = c->E(); // cluster energy
296ea9b4 1644 if(Ec<=0)
1645 return;
6bf90832 1646 Double_t Xc = 0 ; // cluster first moment along X
1647 Double_t Yc = 0 ; // cluster first moment along Y
1648 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
1649 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
1650 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 1651
1652 AliVCaloCells *cells = fEsdCells;
1653 if (!cells)
1654 cells = fAodCells;
1655
6eb6260e 1656 if (!cells)
1657 return;
1658
6bf90832 1659 Int_t ncells = c->GetNCells();
6eb6260e 1660 if (ncells==1)
1661 return;
1662
6eb6260e 1663 for(Int_t j=0; j<ncells; ++j) {
6bf90832 1664 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 1665 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 1666 TVector3 pos;
1667 fGeom->GetGlobal(id,pos);
6eb6260e 1668 Xc += cellen*pos.X();
1669 Yc += cellen*pos.Y();
1670 Sxx += cellen*pos.X()*pos.X();
1671 Syy += cellen*pos.Y()*pos.Y();
1672 Sxy += cellen*pos.X()*pos.Y();
1673 }
1674 Xc /= Ec;
1675 Yc /= Ec;
1676 Sxx /= Ec;
1677 Syy /= Ec;
1678 Sxy /= Ec;
1679 Sxx -= Xc*Xc;
1680 Syy -= Yc*Yc;
1681 Sxy -= Xc*Yc;
f0897c18 1682 Sxx = TMath::Abs(Sxx);
1683 Syy = TMath::Abs(Syy);
296ea9b4 1684 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1685 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1686 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1687 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 1688}
90d5b88b 1689
6bf90832 1690//________________________________________________________________________
0fbe8d4f 1691Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 1692{
1693 // Calculate number of attached cells above emin.
1694
6bf90832 1695 AliVCaloCells *cells = fEsdCells;
1696 if (!cells)
1697 cells = fAodCells;
6bf90832 1698 if (!cells)
296ea9b4 1699 return 0;
6bf90832 1700
296ea9b4 1701 Int_t n = 0;
6bf90832 1702 Int_t ncells = c->GetNCells();
1703 for(Int_t j=0; j<ncells; ++j) {
1704 Int_t id = TMath::Abs(c->GetCellAbsId(j));
1705 Double_t cellen = cells->GetCellAmplitude(id);
1706 if (cellen>=emin)
1707 ++n;
1708 }
1709 return n;
1710}
1711
296ea9b4 1712//________________________________________________________________________
b6c599fe 1713Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 1714{
1715 // Compute isolation based on tracks.
1716
1717 Double_t trkIsolation = 0;
1718 Double_t rad2 = radius*radius;
b6c599fe 1719 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 1720 for(Int_t j = 0; j<ntrks; ++j) {
1721 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1722 if (!track)
1723 continue;
b6c599fe 1724 if (track->Pt()<pt)
1725 continue;
296ea9b4 1726 Float_t eta = track->Eta();
1727 Float_t phi = track->Phi();
0fbe8d4f 1728 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1729 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1730 if(dist>rad2)
1731 continue;
1732 trkIsolation += track->Pt();
1733 }
1734 return trkIsolation;
1735}
3a952328 1736
0fbe8d4f 1737//________________________________________________________________________
1738Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
1739{
1740 // Returns if cluster shared across super modules.
1741
1742 AliVCaloCells *cells = fEsdCells;
1743 if (!cells)
1744 cells = fAodCells;
1745 if (!cells)
1746 return 0;
1747
1748 Int_t n = -1;
1749 Int_t ncells = c->GetNCells();
1750 for(Int_t j=0; j<ncells; ++j) {
1751 Int_t id = TMath::Abs(c->GetCellAbsId(j));
1752 Int_t got = id / (24*48);
1753 if (n==-1) {
1754 n = got;
1755 continue;
1756 }
1757 if (got!=n)
1758 return 1;
1759 }
1760 return 0;
1761}
3a952328 1762