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