]>
Commit | Line | Data |
---|---|---|
ea3fd2d5 | 1 | // $Id$ |
2 | ||
ea3fd2d5 | 3 | #include "AliAnalysisTaskEMCALPi0PbPb.h" |
fa443410 | 4 | #include <TAxis.h> |
717fe7de | 5 | #include <TChain.h> |
6 | #include <TClonesArray.h> | |
f5d4ab70 | 7 | #include <TFile.h> |
296ea9b4 | 8 | #include <TGeoGlobalMagField.h> |
717fe7de | 9 | #include <TH1F.h> |
10 | #include <TH2F.h> | |
11 | #include <TList.h> | |
12 | #include <TLorentzVector.h> | |
f5d4ab70 | 13 | #include <TNtuple.h> |
296ea9b4 | 14 | #include <TProfile.h> |
b3ee6797 | 15 | #include <TString.h> |
296ea9b4 | 16 | #include <TVector2.h> |
717fe7de | 17 | #include "AliAODEvent.h" |
18 | #include "AliAODVertex.h" | |
19 | #include "AliAnalysisManager.h" | |
20 | #include "AliAnalysisTaskEMCALClusterizeFast.h" | |
296ea9b4 | 21 | #include "AliCDBManager.h" |
ea3fd2d5 | 22 | #include "AliCentrality.h" |
ea3fd2d5 | 23 | #include "AliEMCALGeoUtils.h" |
0ec74551 | 24 | #include "AliEMCALGeometry.h" |
296ea9b4 | 25 | #include "AliEMCALRecoUtils.h" |
717fe7de | 26 | #include "AliESDEvent.h" |
27 | #include "AliESDVertex.h" | |
296ea9b4 | 28 | #include "AliESDtrack.h" |
0ec74551 | 29 | #include "AliESDtrackCuts.h" |
296ea9b4 | 30 | #include "AliGeomManager.h" |
b3ee6797 | 31 | #include "AliInputEventHandler.h" |
43cfaa06 | 32 | #include "AliLog.h" |
296ea9b4 | 33 | #include "AliMagF.h" |
0ec74551 | 34 | #include "AliTrackerBase.h" |
ea3fd2d5 | 35 | |
36 | ClassImp(AliAnalysisTaskEMCALPi0PbPb) | |
717fe7de | 37 | |
ea3fd2d5 | 38 | //________________________________________________________________________ |
39 | AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name) | |
40 | : AliAnalysisTaskSE(name), | |
717fe7de | 41 | fCentVar("V0M"), |
42 | fCentFrom(0), | |
43 | fCentTo(100), | |
76332037 | 44 | fVtxZMin(-10), |
45 | fVtxZMax(+10), | |
43cfaa06 | 46 | fUseQualFlag(1), |
717fe7de | 47 | fClusName(), |
f5d4ab70 | 48 | fDoNtuple(0), |
a49742b5 | 49 | fDoAfterburner(0), |
f224d35b | 50 | fAsymMax(1), |
a49742b5 | 51 | fNminCells(1), |
3c661da5 | 52 | fMinE(0.100), |
f224d35b | 53 | fMinErat(0), |
54 | fMinEcc(0), | |
6bf90832 | 55 | fGeoName("EMCAL_FIRSTYEARV1"), |
296ea9b4 | 56 | fMinNClustPerTrack(50), |
0ec74551 | 57 | fMinPtPerTrack(1.0), |
296ea9b4 | 58 | fIsoDist(0.2), |
b3ee6797 | 59 | fTrClassNames(""), |
0ec74551 | 60 | fTrCuts(0), |
2e4d8148 | 61 | fDoTrackMatWithGeom(0), |
62 | fDoConstrain(0), | |
d9f26424 | 63 | fNEvs(0), |
d595acbb | 64 | fGeom(0), |
296ea9b4 | 65 | fReco(0), |
717fe7de | 66 | fOutput(0), |
b3ee6797 | 67 | fTrClassNamesArr(0), |
717fe7de | 68 | fEsdEv(0), |
69 | fAodEv(0), | |
43cfaa06 | 70 | fRecPoints(0), |
717fe7de | 71 | fEsdClusters(0), |
72 | fEsdCells(0), | |
73 | fAodClusters(0), | |
286b47a5 | 74 | fAodCells(0), |
fa443410 | 75 | fPtRanges(0), |
f5d4ab70 | 76 | fNtuple(0), |
296ea9b4 | 77 | fSelTracks(0), |
fa443410 | 78 | fHCuts(0x0), |
79 | fHVertexZ(0x0), | |
76332037 | 80 | fHVertexZ2(0x0), |
fa443410 | 81 | fHCent(0x0), |
82 | fHCentQual(0x0), | |
b3ee6797 | 83 | fHTclsBeforeCuts(0x0), |
84 | fHTclsAfterCuts(0x0), | |
d595acbb | 85 | fHColuRow(0x0), |
86 | fHColuRowE(0x0), | |
87 | fHCellMult(0x0), | |
88 | fHCellE(0x0), | |
89 | fHCellH(0x0), | |
6eb6260e | 90 | fHCellM(0x0), |
a49742b5 | 91 | fHCellM2(0x0), |
296ea9b4 | 92 | fHCellFreqNoCut(0x0), |
2e4d8148 | 93 | fHCellFreqCut100M(0x0), |
94 | fHCellFreqCut300M(0x0), | |
95 | fHCellFreqE(0x0), | |
296ea9b4 | 96 | fHCellCheckE(0x0), |
6eb6260e | 97 | fHClustEccentricity(0), |
fa443410 | 98 | fHClustEtaPhi(0x0), |
99 | fHClustEnergyPt(0x0), | |
100 | fHClustEnergySigma(0x0), | |
d595acbb | 101 | fHClustSigmaSigma(0x0), |
6eb6260e | 102 | fHClustNCellEnergyRatio(0x0), |
fa443410 | 103 | fHPionEtaPhi(0x0), |
104 | fHPionMggPt(0x0), | |
6eb6260e | 105 | fHPionMggAsym(0x0), |
106 | fHPionMggDgg(0x0) | |
ea3fd2d5 | 107 | { |
717fe7de | 108 | // Constructor. |
ea3fd2d5 | 109 | |
d595acbb | 110 | if (!name) |
111 | return; | |
112 | SetName(name); | |
ea3fd2d5 | 113 | DefineInput(0, TChain::Class()); |
ea3fd2d5 | 114 | DefineOutput(1, TList::Class()); |
296ea9b4 | 115 | fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks " |
116 | "AOD:header,vertices,emcalCells,tracks"; | |
ea3fd2d5 | 117 | } |
ea3fd2d5 | 118 | |
717fe7de | 119 | //________________________________________________________________________ |
120 | AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb() | |
121 | { | |
122 | // Destructor. | |
ea3fd2d5 | 123 | |
b3ee6797 | 124 | if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { |
125 | delete fOutput; fOutput = 0; | |
126 | } | |
fa443410 | 127 | delete fPtRanges; fPtRanges = 0; |
d595acbb | 128 | delete fGeom; fGeom = 0; |
296ea9b4 | 129 | delete fReco; fReco = 0; |
b3ee6797 | 130 | delete fTrClassNamesArr; |
296ea9b4 | 131 | delete fSelTracks; |
d595acbb | 132 | delete [] fHColuRow; |
133 | delete [] fHColuRowE; | |
134 | delete [] fHCellMult; | |
296ea9b4 | 135 | delete [] fHCellFreqNoCut; |
2e4d8148 | 136 | delete [] fHCellFreqCut100M; |
137 | delete [] fHCellFreqCut300M; | |
ea3fd2d5 | 138 | } |
717fe7de | 139 | |
ea3fd2d5 | 140 | //________________________________________________________________________ |
141 | void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects() | |
142 | { | |
717fe7de | 143 | // Create user objects here. |
ea3fd2d5 | 144 | |
296ea9b4 | 145 | fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL"); |
146 | fReco = new AliEMCALRecoUtils(); | |
b3ee6797 | 147 | fTrClassNamesArr = fTrClassNames.Tokenize(" "); |
717fe7de | 148 | fOutput = new TList(); |
149 | fOutput->SetOwner(); | |
296ea9b4 | 150 | fSelTracks = new TObjArray; |
ea3fd2d5 | 151 | |
f5d4ab70 | 152 | if (fDoNtuple) { |
153 | TFile *f = OpenFile(1); | |
6bf90832 | 154 | if (f) { |
155 | f->SetCompressionLevel(2); | |
6eb6260e | 156 | fNtuple = new TNtuple(Form("nt%.0fto%.0f",fCentFrom,fCentTo),"nt", |
b3ee6797 | 157 | "run:evt:l0:tcls:cent:pt:eta:phi:e:emax:n:n1:idmax:nsm:db:disp:mn:ms:ecc:sig:tkdz:tkdr:tkep:tkiso:ceiso"); |
6bf90832 | 158 | fNtuple->SetDirectory(f); |
159 | fNtuple->SetAutoFlush(-1024*1024*1024); | |
160 | fNtuple->SetAutoSave(-1024*1024*1024); | |
161 | } | |
f5d4ab70 | 162 | } |
163 | ||
002dcebe | 164 | AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry(); |
165 | Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25; | |
166 | Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25; | |
167 | ||
d595acbb | 168 | // histograms |
a49742b5 | 169 | TH1::SetDefaultSumw2(kTRUE); |
170 | TH2::SetDefaultSumw2(kTRUE); | |
b3ee6797 | 171 | fHCuts = new TH1F("hEventCuts","",5,0.5,5.5); |
172 | fHCuts->GetXaxis()->SetBinLabel(1,"All"); | |
173 | fHCuts->GetXaxis()->SetBinLabel(2,"PS"); | |
174 | fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo)); | |
175 | fHCuts->GetXaxis()->SetBinLabel(4,"QFlag"); | |
176 | fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax)); | |
fa443410 | 177 | fOutput->Add(fHCuts); |
178 | fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25); | |
179 | fHVertexZ->SetXTitle("z [cm]"); | |
180 | fOutput->Add(fHVertexZ); | |
76332037 | 181 | fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25); |
182 | fHVertexZ2->SetXTitle("z [cm]"); | |
183 | fOutput->Add(fHVertexZ2); | |
fa443410 | 184 | fHCent = new TH1F("hCentBeforeCut","",101,-1,100); |
185 | fHCent->SetXTitle(fCentVar.Data()); | |
76332037 | 186 | fOutput->Add(fHCent); |
fa443410 | 187 | fHCentQual = new TH1F("hCentAfterCut","",101,-1,100); |
188 | fHCentQual->SetXTitle(fCentVar.Data()); | |
189 | fOutput->Add(fHCentQual); | |
2e4d8148 | 190 | fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries()); |
191 | fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries()); | |
b3ee6797 | 192 | for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) { |
193 | const char *name = fTrClassNamesArr->At(i)->GetName(); | |
194 | fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name); | |
195 | fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name); | |
196 | } | |
197 | fOutput->Add(fHTclsBeforeCuts); | |
198 | fOutput->Add(fHTclsAfterCuts); | |
90d5b88b | 199 | |
d595acbb | 200 | // histograms for cells |
296ea9b4 | 201 | Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); |
202 | fHColuRow = new TH2*[nsm]; | |
203 | fHColuRowE = new TH2*[nsm]; | |
204 | fHCellMult = new TH1*[nsm]; | |
205 | for (Int_t i = 0; i<nsm; ++i) { | |
2e4d8148 | 206 | fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24); |
d595acbb | 207 | fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i)); |
208 | fHColuRow[i]->SetXTitle("col (i#eta)"); | |
209 | fHColuRow[i]->SetYTitle("row (i#phi)"); | |
2e4d8148 | 210 | fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24); |
90d5b88b | 211 | fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i)); |
d595acbb | 212 | fHColuRowE[i]->SetXTitle("col (i#eta)"); |
213 | fHColuRowE[i]->SetYTitle("row (i#phi)"); | |
214 | fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000); | |
215 | fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i)); | |
90d5b88b | 216 | fHCellMult[i]->SetXTitle("# of cells"); |
d595acbb | 217 | fOutput->Add(fHColuRow[i]); |
218 | fOutput->Add(fHColuRowE[i]); | |
219 | fOutput->Add(fHCellMult[i]); | |
220 | } | |
2e4d8148 | 221 | fHCellE = new TH1F("hCellE","",250,0.,25.); |
d595acbb | 222 | fHCellE->SetXTitle("E_{cell} [GeV]"); |
223 | fOutput->Add(fHCellE); | |
2e4d8148 | 224 | fHCellH = new TH1F ("fHCellHighestE","",250,0.,25.); |
6eb6260e | 225 | fHCellH->SetXTitle("E^{max}_{cell} [GeV]"); |
d595acbb | 226 | fOutput->Add(fHCellH); |
a49742b5 | 227 | fHCellM = new TH1F ("fHCellMeanEperHitCell","",250,0.,2.5); |
6eb6260e | 228 | fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]"); |
229 | fOutput->Add(fHCellM); | |
a49742b5 | 230 | fHCellM2 = new TH1F ("fHCellMeanEperAllCells","",250,0.,1); |
f4ec884e | 231 | fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]"); |
a49742b5 | 232 | fOutput->Add(fHCellM2); |
90d5b88b | 233 | |
2e4d8148 | 234 | fHCellFreqNoCut = new TH1*[nsm]; |
235 | fHCellFreqCut100M = new TH1*[nsm]; | |
236 | fHCellFreqCut300M = new TH1*[nsm]; | |
237 | fHCellFreqE = new TH1*[nsm]; | |
296ea9b4 | 238 | for (Int_t i = 0; i<nsm; ++i){ |
239 | Double_t lbin = i*24*48-0.5; | |
240 | Double_t hbin = lbin+24*48; | |
2e4d8148 | 241 | fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i), |
242 | Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin); | |
243 | fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i), | |
244 | Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin); | |
245 | fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i), | |
246 | Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin); | |
247 | fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i), | |
248 | Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin); | |
296ea9b4 | 249 | fOutput->Add(fHCellFreqNoCut[i]); |
2e4d8148 | 250 | fOutput->Add(fHCellFreqCut100M[i]); |
251 | fOutput->Add(fHCellFreqCut300M[i]); | |
252 | fOutput->Add(fHCellFreqE[i]); | |
296ea9b4 | 253 | } |
254 | if (1) { | |
255 | fHCellCheckE = new TH1*[24*48*nsm]; | |
408dc04c | 256 | memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*)); |
296ea9b4 | 257 | Int_t tcs[1] = {4102}; |
258 | for (UInt_t i = 0; i<sizeof(tcs)/sizeof(Int_t); ++i){ | |
259 | Int_t c=tcs[i]; | |
260 | if (c<24*48*nsm) { | |
261 | fHCellCheckE[i] = new TH1F(Form("fHCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 500, 0, 8); | |
262 | fOutput->Add(fHCellCheckE[i]); | |
263 | } | |
264 | } | |
265 | } | |
266 | ||
d595acbb | 267 | // histograms for clusters |
6eb6260e | 268 | fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1); |
269 | fHClustEccentricity->SetXTitle("#epsilon_{C}"); | |
270 | fOutput->Add(fHClustEccentricity); | |
002dcebe | 271 | fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax); |
fa443410 | 272 | fHClustEtaPhi->SetXTitle("#eta"); |
273 | fHClustEtaPhi->SetYTitle("#varphi"); | |
274 | fOutput->Add(fHClustEtaPhi); | |
275 | fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50); | |
276 | fHClustEnergyPt->SetXTitle("E [GeV]"); | |
6eb6260e | 277 | fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]"); |
fa443410 | 278 | fOutput->Add(fHClustEnergyPt); |
a49742b5 | 279 | fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50); |
280 | fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]"); | |
6eb6260e | 281 | fHClustEnergySigma->SetYTitle("E_{C} [GeV]"); |
d595acbb | 282 | fOutput->Add(fHClustEnergySigma); |
a49742b5 | 283 | fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50); |
6eb6260e | 284 | fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]"); |
285 | fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]"); | |
d595acbb | 286 | fOutput->Add(fHClustSigmaSigma); |
6eb6260e | 287 | fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05); |
288 | fHClustNCellEnergyRatio->SetXTitle("N_{cells}"); | |
289 | fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}"); | |
290 | fOutput->Add(fHClustNCellEnergyRatio); | |
90d5b88b | 291 | |
d595acbb | 292 | // histograms for pion candidates |
002dcebe | 293 | fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax); |
a49742b5 | 294 | fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}"); |
295 | fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}"); | |
fa443410 | 296 | fOutput->Add(fHPionEtaPhi); |
a49742b5 | 297 | fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0); |
fa443410 | 298 | fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
299 | fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]"); | |
300 | fOutput->Add(fHPionMggPt); | |
a49742b5 | 301 | fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1); |
fa443410 | 302 | fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
303 | fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]"); | |
304 | fOutput->Add(fHPionMggAsym); | |
a49742b5 | 305 | fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10); |
6eb6260e | 306 | fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
307 | fHPionMggDgg->SetYTitle("opening angle [grad]"); | |
308 | fOutput->Add(fHPionMggDgg); | |
78204d3d | 309 | const Int_t nbins = 20; |
310 | 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}; | |
fa443410 | 311 | fPtRanges = new TAxis(nbins-1,xbins); |
312 | for (Int_t i = 0; i<=nbins; ++i) { | |
a49742b5 | 313 | fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2); |
fa443410 | 314 | fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
315 | if (i==0) | |
316 | fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0])); | |
317 | else if (i==nbins) | |
318 | fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50")); | |
319 | else | |
320 | fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i])); | |
fa443410 | 321 | fOutput->Add(fHPionInvMasses[i]); |
322 | } | |
296ea9b4 | 323 | |
717fe7de | 324 | PostData(1, fOutput); |
ea3fd2d5 | 325 | } |
326 | ||
327 | //________________________________________________________________________ | |
328 | void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *) | |
329 | { | |
717fe7de | 330 | // Called for each event. |
331 | ||
332 | if (!InputEvent()) | |
333 | return; | |
334 | ||
43cfaa06 | 335 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); |
336 | fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent()); | |
337 | if (fEsdEv) { | |
338 | am->LoadBranch("AliESDRun."); | |
339 | am->LoadBranch("AliESDHeader."); | |
340 | } else { | |
341 | fAodEv = dynamic_cast<AliAODEvent*>(InputEvent()); | |
408dc04c | 342 | if (!fAodEv) { |
343 | AliFatal("Neither ESD nor AOD event found"); | |
344 | return; | |
345 | } | |
43cfaa06 | 346 | am->LoadBranch("header"); |
347 | } | |
348 | ||
b3ee6797 | 349 | Int_t cut = 1; |
350 | fHCuts->Fill(cut++); | |
351 | ||
352 | TString trgclasses; | |
353 | AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader()); | |
354 | if (h) { | |
355 | trgclasses = fEsdEv->GetFiredTriggerClasses(); | |
356 | } else { | |
357 | AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader()); | |
358 | if (h2) | |
359 | trgclasses = h2->GetFiredTriggerClasses(); | |
360 | } | |
361 | for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) { | |
362 | const char *name = fTrClassNamesArr->At(i)->GetName(); | |
363 | if (trgclasses.Contains(name)) | |
364 | fHTclsBeforeCuts->Fill(1+i); | |
365 | } | |
366 | ||
367 | UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); | |
368 | if (res==0) | |
369 | return; | |
370 | ||
371 | if (fHCuts->GetBinContent(2)==0) { | |
2e4d8148 | 372 | if (fDoTrackMatWithGeom && !AliGeomManager::GetGeometry()) { // get geometry |
296ea9b4 | 373 | AliWarning("Accessing geometry from OCDB, this is not very efficient!"); |
374 | AliCDBManager *cdb = AliCDBManager::Instance(); | |
375 | if (!cdb->IsDefaultStorageSet()) | |
376 | cdb->SetDefaultStorage("raw://"); | |
377 | Int_t runno = InputEvent()->GetRunNumber(); | |
378 | if (runno != cdb->GetRun()) | |
379 | cdb->SetRun(runno); | |
380 | AliGeomManager::LoadGeometry(); | |
381 | } | |
382 | ||
002dcebe | 383 | if (!fGeom->GetMatrixForSuperModule(0)) { // set misalignment matrices (stored in first event) |
384 | if (fEsdEv) { | |
385 | for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i) | |
386 | fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i); | |
387 | } else { | |
388 | for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i) | |
389 | fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i); | |
390 | } | |
d595acbb | 391 | } |
296ea9b4 | 392 | |
393 | if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map | |
394 | if (fEsdEv) { | |
395 | const AliESDRun *erun = fEsdEv->GetESDRun(); | |
396 | AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(), | |
397 | erun->GetCurrentDip(), | |
398 | AliMagF::kConvLHC, | |
399 | kFALSE, | |
400 | erun->GetBeamEnergy(), | |
401 | erun->GetBeamType()); | |
402 | TGeoGlobalMagField::Instance()->SetField(field); | |
403 | } else { | |
404 | Double_t pol = -1; //polarity | |
405 | Double_t be = -1; //beam energy | |
406 | AliMagF::BeamType_t btype = AliMagF::kBeamTypepp; | |
407 | Int_t runno = fAodEv->GetRunNumber(); | |
408 | if (runno>=136851 && runno<138275) { | |
409 | pol = -1; | |
410 | be = 2760; | |
411 | btype = AliMagF::kBeamTypeAA; | |
412 | } else if (runno>=138275 && runno<=139517) { | |
413 | pol = +1; | |
414 | be = 2760; | |
415 | btype = AliMagF::kBeamTypeAA; | |
416 | } else { | |
417 | AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno)); | |
418 | } | |
419 | TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be)); | |
420 | } | |
421 | } | |
d595acbb | 422 | } |
423 | ||
fa443410 | 424 | fHCuts->Fill(cut++); |
286b47a5 | 425 | |
717fe7de | 426 | const AliCentrality *centP = InputEvent()->GetCentrality(); |
427 | Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar); | |
fa443410 | 428 | fHCent->Fill(cent); |
717fe7de | 429 | if (cent<fCentFrom||cent>fCentTo) |
286b47a5 | 430 | return; |
43cfaa06 | 431 | |
fa443410 | 432 | fHCuts->Fill(cut++); |
286b47a5 | 433 | |
43cfaa06 | 434 | if (fUseQualFlag) { |
435 | if (centP->GetQuality()>0) | |
436 | return; | |
437 | } | |
286b47a5 | 438 | |
fa443410 | 439 | fHCentQual->Fill(cent); |
440 | fHCuts->Fill(cut++); | |
717fe7de | 441 | |
43cfaa06 | 442 | if (fEsdEv) { |
fa443410 | 443 | am->LoadBranch("PrimaryVertex."); |
444 | am->LoadBranch("SPDVertex."); | |
445 | am->LoadBranch("TPCVertex."); | |
43cfaa06 | 446 | } else { |
447 | fAodEv = dynamic_cast<AliAODEvent*>(InputEvent()); | |
448 | am->LoadBranch("vertices"); | |
3c661da5 | 449 | if (!fAodEv) return; |
43cfaa06 | 450 | } |
451 | ||
452 | const AliVVertex *vertex = InputEvent()->GetPrimaryVertex(); | |
717fe7de | 453 | if (!vertex) |
454 | return; | |
455 | ||
fa443410 | 456 | fHVertexZ->Fill(vertex->GetZ()); |
286b47a5 | 457 | |
717fe7de | 458 | if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax) |
459 | return; | |
286b47a5 | 460 | |
fa443410 | 461 | fHCuts->Fill(cut++); |
76332037 | 462 | fHVertexZ2->Fill(vertex->GetZ()); |
717fe7de | 463 | |
b3ee6797 | 464 | // count number of accepted events |
465 | ++fNEvs; | |
466 | ||
467 | for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) { | |
468 | const char *name = fTrClassNamesArr->At(i)->GetName(); | |
469 | if (trgclasses.Contains(name)) | |
470 | fHTclsAfterCuts->Fill(1+i); | |
471 | } | |
472 | ||
43cfaa06 | 473 | fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used |
474 | fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached | |
475 | fEsdCells = 0; // will be set if ESD input used | |
476 | fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached | |
477 | // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode | |
478 | fAodCells = 0; // will be set if AOD input used | |
717fe7de | 479 | |
43cfaa06 | 480 | // deal with special output from AliAnalysisTaskEMCALClusterizeFast first |
481 | Bool_t clusattached = 0; | |
482 | Bool_t recalibrated = 0; | |
483 | if (1 && !fClusName.IsNull()) { | |
484 | AliAnalysisTaskEMCALClusterizeFast *cltask = 0; | |
485 | TObjArray *ts = am->GetTasks(); | |
486 | cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName)); | |
487 | if (cltask && cltask->GetClusters()) { | |
d595acbb | 488 | fRecPoints = const_cast<TObjArray*>(cltask->GetClusters()); |
43cfaa06 | 489 | clusattached = cltask->GetAttachClusters(); |
490 | if (cltask->GetCalibData()!=0) | |
491 | recalibrated = kTRUE; | |
492 | } | |
493 | } | |
494 | if (1 && AODEvent() && !fClusName.IsNull()) { | |
717fe7de | 495 | TList *l = AODEvent()->GetList(); |
496 | TClonesArray *clus = 0; | |
497 | if (l) { | |
498 | clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName)); | |
499 | fAodClusters = clus; | |
500 | } | |
ea3fd2d5 | 501 | } |
717fe7de | 502 | |
503 | if (fEsdEv) { // ESD input mode | |
43cfaa06 | 504 | if (1 && (!fRecPoints||clusattached)) { |
505 | if (!clusattached) | |
506 | am->LoadBranch("CaloClusters"); | |
717fe7de | 507 | TList *l = fEsdEv->GetList(); |
508 | TClonesArray *clus = 0; | |
509 | if (l) { | |
510 | clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters")); | |
511 | fEsdClusters = clus; | |
ea3fd2d5 | 512 | } |
513 | } | |
43cfaa06 | 514 | if (1) { |
515 | if (!recalibrated) | |
516 | am->LoadBranch("EMCALCells."); | |
717fe7de | 517 | fEsdCells = fEsdEv->GetEMCALCells(); |
518 | } | |
519 | } else if (fAodEv) { // AOD input mode | |
43cfaa06 | 520 | if (1 && (!fAodClusters || clusattached)) { |
521 | if (!clusattached) | |
522 | am->LoadBranch("caloClusters"); | |
717fe7de | 523 | TList *l = fAodEv->GetList(); |
524 | TClonesArray *clus = 0; | |
525 | if (l) { | |
526 | clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters")); | |
527 | fAodClusters = clus; | |
ea3fd2d5 | 528 | } |
529 | } | |
43cfaa06 | 530 | if (1) { |
531 | if (!recalibrated) | |
532 | am->LoadBranch("emcalCells"); | |
717fe7de | 533 | fAodCells = fAodEv->GetEMCALCells(); |
534 | } | |
535 | } else { | |
536 | AliFatal("Impossible to not have either pointer to ESD or AOD event"); | |
537 | } | |
ea3fd2d5 | 538 | |
43cfaa06 | 539 | if (1) { |
540 | AliDebug(2,Form("fRecPoints set: %p", fRecPoints)); | |
541 | AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters)); | |
542 | AliDebug(2,Form("fEsdCells set: %p", fEsdCells)); | |
543 | AliDebug(2,Form("fAodClusters set: %p", fAodClusters)); | |
544 | AliDebug(2,Form("fAodCells set: %p", fAodCells)); | |
545 | } | |
546 | ||
a49742b5 | 547 | if (fDoAfterburner) |
548 | ClusterAfterburner(); | |
6eb6260e | 549 | |
296ea9b4 | 550 | CalcTracks(); |
551 | CalcClusterProps(); | |
552 | ||
286b47a5 | 553 | FillCellHists(); |
554 | FillClusHists(); | |
555 | FillPionHists(); | |
323834f0 | 556 | FillOtherHists(); |
ea3fd2d5 | 557 | |
717fe7de | 558 | PostData(1, fOutput); |
ea3fd2d5 | 559 | } |
560 | ||
561 | //________________________________________________________________________ | |
562 | void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *) | |
563 | { | |
717fe7de | 564 | // Terminate called at the end of analysis. |
f5d4ab70 | 565 | |
566 | if (fNtuple) { | |
567 | TFile *f = OpenFile(1); | |
568 | if (f) | |
569 | fNtuple->Write(); | |
570 | } | |
d9f26424 | 571 | |
b3ee6797 | 572 | AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs)); |
ea3fd2d5 | 573 | } |
ea3fd2d5 | 574 | |
296ea9b4 | 575 | //________________________________________________________________________ |
576 | void AliAnalysisTaskEMCALPi0PbPb::CalcTracks() | |
577 | { | |
578 | // Calculate track properties. | |
579 | ||
580 | fSelTracks->Clear(); | |
581 | ||
582 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); | |
583 | TClonesArray *tracks = 0; | |
584 | if (fEsdEv) { | |
585 | am->LoadBranch("Tracks"); | |
586 | TList *l = fEsdEv->GetList(); | |
587 | tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks")); | |
588 | } else { | |
589 | am->LoadBranch("tracks"); | |
590 | TList *l = fAodEv->GetList(); | |
591 | tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks")); | |
592 | } | |
593 | ||
594 | if (!tracks) | |
595 | return; | |
596 | ||
0ec74551 | 597 | AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry(); |
598 | Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25; | |
599 | Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25; | |
600 | ||
296ea9b4 | 601 | if (fEsdEv) { |
2e4d8148 | 602 | if (fDoConstrain) |
603 | fSelTracks->SetOwner(kTRUE); | |
0ec74551 | 604 | am->LoadBranch("PrimaryVertex."); |
605 | const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD(); | |
606 | am->LoadBranch("SPDVertex."); | |
607 | const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD(); | |
608 | am->LoadBranch("Tracks"); | |
609 | const Int_t Ntracks = tracks->GetEntries(); | |
610 | for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { | |
611 | AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks)); | |
612 | if (!track) { | |
613 | AliWarning(Form("Could not receive track %d\n", iTracks)); | |
614 | continue; | |
615 | } | |
616 | if (fTrCuts && !fTrCuts->IsSelected(track)) | |
296ea9b4 | 617 | continue; |
618 | Double_t eta = track->Eta(); | |
619 | if (eta<-1||eta>1) | |
620 | continue; | |
0ec74551 | 621 | if (track->Phi()<phimin||track->Phi()>phimax) |
622 | continue; | |
296ea9b4 | 623 | if(track->GetTPCNcls()<fMinNClustPerTrack) |
624 | continue; | |
2e4d8148 | 625 | |
626 | if (!fDoConstrain) { | |
627 | fSelTracks->Add(track); | |
628 | continue; | |
629 | } | |
630 | ||
0ec74551 | 631 | AliESDtrack copyt(*track); |
632 | Double_t bfield[3]; | |
633 | copyt.GetBxByBz(bfield); | |
634 | AliExternalTrackParam tParam; | |
635 | Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam); | |
636 | if (!relate) | |
637 | continue; | |
638 | copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance()); | |
2e4d8148 | 639 | |
0ec74551 | 640 | Double_t p[3] = { 0. }; |
641 | copyt.GetPxPyPz(p); | |
642 | Double_t pos[3] = { 0. }; | |
643 | copyt.GetXYZ(pos); | |
644 | Double_t covTr[21] = { 0. }; | |
645 | copyt.GetCovarianceXYZPxPyPz(covTr); | |
646 | Double_t pid[10] = { 0. }; | |
647 | copyt.GetESDpid(pid); | |
648 | AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(), | |
649 | copyt.GetLabel(), | |
650 | p, | |
651 | kTRUE, | |
652 | pos, | |
653 | kFALSE, | |
654 | covTr, | |
655 | (Short_t)copyt.GetSign(), | |
656 | copyt.GetITSClusterMap(), | |
657 | pid, | |
658 | 0,/*fPrimaryVertex,*/ | |
659 | kTRUE, // check if this is right | |
660 | vtx->UsesTrack(copyt.GetID())); | |
661 | aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap()); | |
662 | aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap()); | |
663 | Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ; | |
664 | if(ndf>0) | |
665 | aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf); | |
666 | else | |
667 | aTrack->SetChi2perNDF(-1); | |
668 | aTrack->SetFlags(copyt.GetStatus()); | |
669 | aTrack->SetTPCPointsF(copyt.GetTPCNclsF()); | |
670 | fSelTracks->Add(aTrack); | |
296ea9b4 | 671 | } |
672 | } else { | |
673 | Int_t ntracks = tracks->GetEntries(); | |
674 | for (Int_t i=0; i<ntracks; ++i) { | |
675 | AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i)); | |
676 | if (!track) | |
677 | continue; | |
678 | Double_t eta = track->Eta(); | |
679 | if (eta<-1||eta>1) | |
680 | continue; | |
0ec74551 | 681 | if (track->Phi()<phimin||track->Phi()>phimax) |
682 | continue; | |
296ea9b4 | 683 | if(track->GetTPCNcls()<fMinNClustPerTrack) |
684 | continue; | |
c4236a58 | 685 | |
71fed3bd | 686 | if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL |
c4236a58 | 687 | AliExternalTrackParam tParam(track); |
688 | if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) { | |
689 | Double_t trkPos[3]; | |
690 | tParam.GetXYZ(trkPos); | |
691 | track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]); | |
692 | } | |
693 | } | |
296ea9b4 | 694 | fSelTracks->Add(track); |
695 | } | |
696 | } | |
697 | } | |
698 | ||
699 | //________________________________________________________________________ | |
700 | void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps() | |
701 | { | |
702 | // Calculate cluster properties | |
703 | ||
704 | TObjArray *clusters = fEsdClusters; | |
705 | if (!clusters) | |
706 | clusters = fAodClusters; | |
707 | if (!clusters) | |
708 | return; | |
709 | ||
710 | Int_t nclus = clusters->GetEntries(); | |
c4236a58 | 711 | Int_t ntrks = fSelTracks->GetEntries(); |
712 | ||
296ea9b4 | 713 | for(Int_t i = 0; i<nclus; ++i) { |
714 | fClusProps[i].Reset(); | |
715 | ||
716 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); | |
717 | if (!clus) | |
718 | continue; | |
719 | if (!clus->IsEMCAL()) | |
720 | continue; | |
721 | if (clus->E()<fMinE) | |
722 | continue; | |
723 | ||
724 | Float_t clsPos[3] = {0}; | |
725 | clus->GetPosition(clsPos); | |
04785f50 | 726 | TVector3 clsVec(clsPos); |
c4236a58 | 727 | Double_t vertex[3] = {0}; |
728 | InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); | |
729 | TLorentzVector clusterVec; | |
730 | clus->GetMomentum(clusterVec,vertex); | |
731 | Double_t clsEta = clusterVec.Eta(); | |
296ea9b4 | 732 | |
296ea9b4 | 733 | Double_t mind2 = 1e10; |
734 | for(Int_t j = 0; j<ntrks; ++j) { | |
735 | AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j)); | |
736 | if (!track) | |
737 | continue; | |
71fed3bd | 738 | Double_t pt = track->Pt(); |
739 | if (pt<fMinPtPerTrack) | |
296ea9b4 | 740 | continue; |
2e4d8148 | 741 | if (TMath::Abs(clsEta-track->Eta())>0.5) |
c4236a58 | 742 | continue; |
743 | ||
296ea9b4 | 744 | Float_t tmpR=-1, tmpZ=-1; |
2e4d8148 | 745 | if (!fDoTrackMatWithGeom) { |
746 | ||
747 | AliExternalTrackParam *tParam = 0; | |
748 | if (fEsdEv) { | |
749 | AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track); | |
750 | tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam()); | |
751 | } else | |
752 | tParam = new AliExternalTrackParam(track); | |
753 | ||
754 | Double_t bfield[3]; | |
755 | track->GetBxByBz(bfield); | |
04785f50 | 756 | TVector3 vec(clsPos); |
757 | Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad(); | |
758 | vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system | |
2e4d8148 | 759 | tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system |
760 | Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield); | |
761 | if (!ret) { | |
762 | delete tParam; | |
04785f50 | 763 | continue; |
2e4d8148 | 764 | } |
04785f50 | 765 | Double_t trkPos[3]; |
2e4d8148 | 766 | tParam->GetXYZ(trkPos); //Get the extrapolated global position |
04785f50 | 767 | tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) ); |
768 | tmpZ = clsPos[2]-trkPos[2]; | |
2e4d8148 | 769 | delete tParam; |
770 | } else { | |
771 | if (TMath::Abs(clsEta-track->Eta())>fIsoDist) | |
772 | continue; | |
773 | AliExternalTrackParam tParam(track); | |
774 | if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ)) | |
775 | continue; | |
04785f50 | 776 | } |
2e4d8148 | 777 | |
c4236a58 | 778 | Double_t d2 = tmpR; |
296ea9b4 | 779 | if (mind2>d2) { |
780 | mind2=d2; | |
781 | fClusProps[i].fTrIndex = j; | |
782 | fClusProps[i].fTrDz = tmpZ; | |
c4236a58 | 783 | fClusProps[i].fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ); |
296ea9b4 | 784 | fClusProps[i].fTrDist = d2; |
785 | fClusProps[i].fTrEp = clus->E()/track->P(); | |
786 | } | |
787 | } | |
788 | ||
c4236a58 | 789 | if (0 && (fClusProps[i].fTrIndex>=0)) { |
790 | cout << i << " " << fClusProps[i].fTrIndex << ": Dr " << fClusProps[i].fTrDr << " " << " Dz " << fClusProps[i].fTrDz << endl; | |
791 | } | |
2e4d8148 | 792 | |
04785f50 | 793 | fClusProps[i].fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist); |
794 | fClusProps[i].fTrLowPtIso = 0; | |
795 | fClusProps[i].fCellIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist); | |
296ea9b4 | 796 | } |
797 | } | |
798 | ||
323834f0 | 799 | //________________________________________________________________________ |
800 | void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner() | |
801 | { | |
296ea9b4 | 802 | // Run custer reconstruction afterburner. |
323834f0 | 803 | |
804 | AliVCaloCells *cells = fEsdCells; | |
805 | if (!cells) | |
806 | cells = fAodCells; | |
807 | ||
808 | if (!cells) | |
809 | return; | |
810 | ||
811 | Int_t ncells = cells->GetNumberOfCells(); | |
812 | if (ncells<=0) | |
813 | return; | |
814 | ||
815 | Double_t cellMeanE = 0, cellSigE = 0; | |
816 | for (Int_t i = 0; i<ncells; ++i) { | |
817 | Double_t cellE = cells->GetAmplitude(i); | |
818 | cellMeanE += cellE; | |
819 | cellSigE += cellE*cellE; | |
820 | } | |
821 | cellMeanE /= ncells; | |
822 | cellSigE /= ncells; | |
823 | cellSigE -= cellMeanE*cellMeanE; | |
824 | if (cellSigE<0) | |
825 | cellSigE = 0; | |
826 | cellSigE = TMath::Sqrt(cellSigE / ncells); | |
827 | ||
828 | Double_t subE = cellMeanE - 7*cellSigE; | |
829 | if (subE<0) | |
830 | return; | |
831 | ||
832 | for (Int_t i = 0; i<ncells; ++i) { | |
833 | Short_t id=-1; | |
834 | Double_t amp=0,time=0; | |
835 | if (!cells->GetCell(i, id, amp, time)) | |
836 | continue; | |
837 | amp -= cellMeanE; | |
838 | if (amp<0.001) | |
839 | amp = 0; | |
840 | cells->SetCell(i, id, amp, time); | |
841 | } | |
842 | ||
843 | TObjArray *clusters = fEsdClusters; | |
844 | if (!clusters) | |
845 | clusters = fAodClusters; | |
323834f0 | 846 | if (!clusters) |
847 | return; | |
848 | ||
849 | Int_t nclus = clusters->GetEntries(); | |
850 | for (Int_t i = 0; i<nclus; ++i) { | |
851 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); | |
852 | if (!clus->IsEMCAL()) | |
853 | continue; | |
854 | Int_t nc = clus->GetNCells(); | |
855 | Double_t clusE = 0; | |
856 | UShort_t ids[100] = {0}; | |
857 | Double_t fra[100] = {0}; | |
858 | for (Int_t j = 0; j<nc; ++j) { | |
859 | Short_t id = TMath::Abs(clus->GetCellAbsId(j)); | |
860 | Double_t cen = cells->GetCellAmplitude(id); | |
861 | clusE += cen; | |
862 | if (cen>0) { | |
863 | ids[nc] = id; | |
864 | ++nc; | |
865 | } | |
866 | } | |
867 | if (clusE<=0) { | |
868 | clusters->RemoveAt(i); | |
869 | continue; | |
870 | } | |
871 | ||
872 | for (Int_t j = 0; j<nc; ++j) { | |
873 | Short_t id = ids[j]; | |
874 | Double_t cen = cells->GetCellAmplitude(id); | |
875 | fra[j] = cen/clusE; | |
876 | } | |
877 | clus->SetE(clusE); | |
878 | AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus); | |
879 | if (aodclus) { | |
880 | aodclus->Clear(""); | |
881 | aodclus->SetNCells(nc); | |
882 | aodclus->SetCellsAmplitudeFraction(fra); | |
883 | aodclus->SetCellsAbsId(ids); | |
884 | } | |
885 | } | |
886 | clusters->Compress(); | |
887 | } | |
888 | ||
286b47a5 | 889 | //________________________________________________________________________ |
890 | void AliAnalysisTaskEMCALPi0PbPb::FillCellHists() | |
891 | { | |
892 | // Fill histograms related to cell properties. | |
d595acbb | 893 | |
90d5b88b | 894 | AliVCaloCells *cells = fEsdCells; |
895 | if (!cells) | |
896 | cells = fAodCells; | |
897 | ||
296ea9b4 | 898 | if (!cells) |
899 | return; | |
90d5b88b | 900 | |
296ea9b4 | 901 | Int_t cellModCount[12] = {0}; |
902 | Double_t cellMaxE = 0; | |
903 | Double_t cellMeanE = 0; | |
904 | Int_t ncells = cells->GetNumberOfCells(); | |
905 | if (ncells==0) | |
906 | return; | |
907 | ||
908 | Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); | |
909 | ||
910 | for (Int_t i = 0; i<ncells; ++i) { | |
911 | Int_t absID = TMath::Abs(cells->GetCellNumber(i)); | |
912 | Double_t cellE = cells->GetAmplitude(i); | |
913 | fHCellE->Fill(cellE); | |
914 | if (cellE>cellMaxE) | |
915 | cellMaxE = cellE; | |
916 | cellMeanE += cellE; | |
917 | ||
918 | Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1; | |
919 | Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta); | |
920 | if (!ret) { | |
921 | AliError(Form("Could not get cell index for %d", absID)); | |
922 | continue; | |
923 | } | |
924 | ++cellModCount[iSM]; | |
925 | Int_t iPhi=-1, iEta=-1; | |
926 | fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta); | |
927 | fHColuRow[iSM]->Fill(iEta,iPhi,1); | |
928 | fHColuRowE[iSM]->Fill(iEta,iPhi,cellE); | |
929 | fHCellFreqNoCut[iSM]->Fill(absID); | |
2e4d8148 | 930 | if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID); |
931 | if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID); | |
296ea9b4 | 932 | if (fHCellCheckE && fHCellCheckE[absID]) |
933 | fHCellCheckE[absID]->Fill(cellE); | |
2e4d8148 | 934 | fHCellFreqE[iSM]->Fill(absID, cellE); |
296ea9b4 | 935 | } |
936 | fHCellH->Fill(cellMaxE); | |
937 | cellMeanE /= ncells; | |
938 | fHCellM->Fill(cellMeanE); | |
939 | fHCellM2->Fill(cellMeanE*ncells/24/48/nsm); | |
940 | for (Int_t i=0; i<nsm; ++i) | |
941 | fHCellMult[i]->Fill(cellModCount[i]); | |
286b47a5 | 942 | } |
943 | ||
944 | //________________________________________________________________________ | |
945 | void AliAnalysisTaskEMCALPi0PbPb::FillClusHists() | |
946 | { | |
90d5b88b | 947 | // Fill histograms related to cluster properties. |
fa443410 | 948 | |
90d5b88b | 949 | TObjArray *clusters = fEsdClusters; |
950 | if (!clusters) | |
951 | clusters = fAodClusters; | |
296ea9b4 | 952 | if (!clusters) |
953 | return; | |
90d5b88b | 954 | |
296ea9b4 | 955 | Double_t vertex[3] = {0}; |
956 | InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); | |
957 | ||
958 | Int_t nclus = clusters->GetEntries(); | |
959 | for(Int_t i = 0; i<nclus; ++i) { | |
960 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); | |
961 | if (!clus) | |
962 | continue; | |
963 | if (!clus->IsEMCAL()) | |
964 | continue; | |
90d5b88b | 965 | TLorentzVector clusterVec; |
296ea9b4 | 966 | clus->GetMomentum(clusterVec,vertex); |
967 | Double_t maxAxis = 0, minAxis = 0; | |
968 | GetSigma(clus,maxAxis,minAxis); | |
969 | Double_t clusterEcc = 0; | |
970 | if (maxAxis > 0) | |
971 | clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis)); | |
972 | clus->SetChi2(clusterEcc); // store ecc in chi2 | |
973 | fHClustEccentricity->Fill(clusterEcc); | |
974 | fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi()); | |
975 | fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt()); | |
976 | fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E()); | |
977 | fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis); | |
978 | fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E()); | |
979 | if (fNtuple) { | |
980 | if (clus->E()<fMinE) | |
fa443410 | 981 | continue; |
b3ee6797 | 982 | Float_t vals[25]; |
002dcebe | 983 | TString trgclasses; |
296ea9b4 | 984 | vals[0] = InputEvent()->GetRunNumber(); |
985 | vals[1] = (((UInt_t)InputEvent()->GetOrbitNumber() << 12) | (UInt_t)InputEvent()->GetBunchCrossNumber()); | |
986 | if (vals[1]<=0) | |
987 | vals[1] = fNEvs; | |
988 | AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader()); | |
002dcebe | 989 | if (h) { |
296ea9b4 | 990 | vals[2] = h->GetL0TriggerInputs(); |
002dcebe | 991 | trgclasses = fEsdEv->GetFiredTriggerClasses(); |
992 | } | |
296ea9b4 | 993 | else { |
994 | AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader()); | |
002dcebe | 995 | if (h2) { |
296ea9b4 | 996 | vals[2] = h2->GetL0TriggerInputs(); |
002dcebe | 997 | trgclasses = h2->GetFiredTriggerClasses(); |
998 | } else | |
296ea9b4 | 999 | vals[2] = 0; |
a49742b5 | 1000 | } |
002dcebe | 1001 | vals[3] = 0; |
b3ee6797 | 1002 | |
1003 | for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) { | |
1004 | const char *name = fTrClassNamesArr->At(j)->GetName(); | |
1005 | if (trgclasses.Contains(name)) | |
1006 | vals[3] += TMath::Power(2,j); | |
1007 | } | |
002dcebe | 1008 | vals[4] = InputEvent()->GetCentrality()->GetCentralityPercentileUnchecked(fCentVar); |
1009 | vals[5] = clusterVec.Pt(); | |
1010 | vals[6] = clusterVec.Eta(); | |
1011 | vals[7] = clusterVec.Phi(); | |
1012 | vals[8] = clusterVec.E(); | |
1013 | vals[9] = GetMaxCellEnergy(clus); | |
b3ee6797 | 1014 | vals[10] = clus->GetNCells(); |
002dcebe | 1015 | vals[11] = GetNCells(clus,0.100); |
b3ee6797 | 1016 | vals[12] = clus->GetCellAbsId(0); |
1017 | vals[13] = fGeom->GetSuperModuleNumber(clus->GetCellAbsId(0)); | |
1018 | vals[14] = clus->GetDistanceToBadChannel(); | |
1019 | vals[15] = clus->GetDispersion(); | |
1020 | vals[16] = clus->GetM20(); | |
1021 | vals[17] = clus->GetM02(); | |
1022 | vals[18] = clusterEcc; | |
1023 | vals[19] = maxAxis; | |
1024 | vals[20] = fClusProps[i].fTrDz; | |
1025 | vals[21] = fClusProps[i].fTrDr; | |
1026 | vals[22] = fClusProps[i].fTrEp; | |
1027 | vals[23] = fClusProps[i].fTrIso; | |
1028 | vals[24] = fClusProps[i].fCellIso; | |
296ea9b4 | 1029 | fNtuple->Fill(vals); |
fa443410 | 1030 | } |
1031 | } | |
286b47a5 | 1032 | } |
1033 | ||
1034 | //________________________________________________________________________ | |
1035 | void AliAnalysisTaskEMCALPi0PbPb::FillPionHists() | |
1036 | { | |
1037 | // Fill histograms related to pions. | |
286b47a5 | 1038 | |
296ea9b4 | 1039 | Double_t vertex[3] = {0}; |
fa443410 | 1040 | InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); |
1041 | ||
90d5b88b | 1042 | TObjArray *clusters = fEsdClusters; |
1043 | if (!clusters) | |
1044 | clusters = fAodClusters; | |
fa443410 | 1045 | |
90d5b88b | 1046 | if (clusters) { |
1047 | TLorentzVector clusterVec1; | |
1048 | TLorentzVector clusterVec2; | |
1049 | TLorentzVector pionVec; | |
1050 | ||
1051 | Int_t nclus = clusters->GetEntries(); | |
fa443410 | 1052 | for (Int_t i = 0; i<nclus; ++i) { |
90d5b88b | 1053 | AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i)); |
fa443410 | 1054 | if (!clus1) |
1055 | continue; | |
1056 | if (!clus1->IsEMCAL()) | |
1057 | continue; | |
3c661da5 | 1058 | if (clus1->E()<fMinE) |
6eb6260e | 1059 | continue; |
a49742b5 | 1060 | if (clus1->GetNCells()<fNminCells) |
1061 | continue; | |
f224d35b | 1062 | if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat) |
1063 | continue; | |
3c661da5 | 1064 | if (clus1->Chi2()<fMinEcc) // eccentricity cut |
f224d35b | 1065 | continue; |
fa443410 | 1066 | clus1->GetMomentum(clusterVec1,vertex); |
1067 | for (Int_t j = i+1; j<nclus; ++j) { | |
90d5b88b | 1068 | AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j)); |
fa443410 | 1069 | if (!clus2) |
1070 | continue; | |
1071 | if (!clus2->IsEMCAL()) | |
1072 | continue; | |
3c661da5 | 1073 | if (clus2->E()<fMinE) |
6eb6260e | 1074 | continue; |
a49742b5 | 1075 | if (clus2->GetNCells()<fNminCells) |
1076 | continue; | |
f224d35b | 1077 | if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat) |
1078 | continue; | |
3c661da5 | 1079 | if (clus2->Chi2()<fMinEcc) // eccentricity cut |
f224d35b | 1080 | continue; |
fa443410 | 1081 | clus2->GetMomentum(clusterVec2,vertex); |
1082 | pionVec = clusterVec1 + clusterVec2; | |
1083 | Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E(); | |
6eb6260e | 1084 | Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect()); |
d595acbb | 1085 | if (pionZgg < fAsymMax) { |
1086 | fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi()); | |
1087 | fHPionMggPt->Fill(pionVec.M(),pionVec.Pt()); | |
1088 | fHPionMggAsym->Fill(pionVec.M(),pionZgg); | |
6eb6260e | 1089 | fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle); |
d595acbb | 1090 | Int_t bin = fPtRanges->FindBin(pionVec.Pt()); |
1091 | fHPionInvMasses[bin]->Fill(pionVec.M()); | |
1092 | } | |
fa443410 | 1093 | } |
1094 | } | |
90d5b88b | 1095 | } |
fa443410 | 1096 | } |
d595acbb | 1097 | |
6eb6260e | 1098 | //________________________________________________________________________ |
323834f0 | 1099 | void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists() |
6eb6260e | 1100 | { |
323834f0 | 1101 | // Fill histograms related to cell properties. |
6eb6260e | 1102 | } |
1103 | ||
d595acbb | 1104 | //________________________________________________________________________ |
296ea9b4 | 1105 | Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const |
1106 | { | |
1107 | // Compute isolation based on cell content. | |
1108 | ||
1109 | AliVCaloCells *cells = fEsdCells; | |
1110 | if (!cells) | |
1111 | cells = fAodCells; | |
1112 | if (!cells) | |
1113 | return 0; | |
1114 | ||
1115 | Double_t cellIsolation = 0; | |
1116 | Double_t rad2 = radius*radius; | |
1117 | Int_t ncells = cells->GetNumberOfCells(); | |
1118 | for (Int_t i = 0; i<ncells; ++i) { | |
1119 | Int_t absID = TMath::Abs(cells->GetCellNumber(i)); | |
1120 | Double_t cellE = cells->GetAmplitude(i); | |
1121 | Float_t eta, phi; | |
1122 | fGeom->EtaPhiFromIndex(absID,eta,phi); | |
1123 | Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi); | |
1124 | Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff; | |
1125 | if(dist>rad2) | |
1126 | continue; | |
1127 | cellIsolation += cellE; | |
1128 | } | |
1129 | return cellIsolation; | |
1130 | } | |
1131 | ||
1132 | //________________________________________________________________________ | |
1133 | Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster) const | |
d595acbb | 1134 | { |
90d5b88b | 1135 | // Get maximum energy of attached cell. |
1136 | ||
1137 | Double_t maxe = 0; | |
90d5b88b | 1138 | Int_t ncells = cluster->GetNCells(); |
1139 | if (fEsdCells) { | |
1140 | for (Int_t i=0; i<ncells; i++) { | |
27422847 | 1141 | Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); |
f0897c18 | 1142 | if (e>maxe) { |
90d5b88b | 1143 | maxe = e; |
f0897c18 | 1144 | } |
90d5b88b | 1145 | } |
1146 | } else { | |
1147 | for (Int_t i=0; i<ncells; i++) { | |
27422847 | 1148 | Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); |
90d5b88b | 1149 | if (e>maxe) |
1150 | maxe = e; | |
1151 | } | |
1152 | } | |
6eb6260e | 1153 | return maxe; |
d595acbb | 1154 | } |
1155 | ||
1156 | //________________________________________________________________________ | |
296ea9b4 | 1157 | void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const |
d595acbb | 1158 | { |
90d5b88b | 1159 | // Calculate the (E) weighted variance along the longer (eigen) axis. |
1160 | ||
6bf90832 | 1161 | sigmaMax = 0; // cluster variance along its longer axis |
1162 | sigmaMin = 0; // cluster variance along its shorter axis | |
1163 | Double_t Ec = c->E(); // cluster energy | |
296ea9b4 | 1164 | if(Ec<=0) |
1165 | return; | |
6bf90832 | 1166 | Double_t Xc = 0 ; // cluster first moment along X |
1167 | Double_t Yc = 0 ; // cluster first moment along Y | |
1168 | Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2) | |
1169 | Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2) | |
1170 | Double_t Syy = 0 ; // cluster covariance^2 | |
90d5b88b | 1171 | |
1172 | AliVCaloCells *cells = fEsdCells; | |
1173 | if (!cells) | |
1174 | cells = fAodCells; | |
1175 | ||
6eb6260e | 1176 | if (!cells) |
1177 | return; | |
1178 | ||
6bf90832 | 1179 | Int_t ncells = c->GetNCells(); |
6eb6260e | 1180 | if (ncells==1) |
1181 | return; | |
1182 | ||
1183 | TVector3 pos; | |
1184 | for(Int_t j=0; j<ncells; ++j) { | |
6bf90832 | 1185 | Int_t id = TMath::Abs(c->GetCellAbsId(j)); |
27422847 | 1186 | fGeom->GetGlobal(id,pos); |
6eb6260e | 1187 | Double_t cellen = cells->GetCellAmplitude(id); |
1188 | Xc += cellen*pos.X(); | |
1189 | Yc += cellen*pos.Y(); | |
1190 | Sxx += cellen*pos.X()*pos.X(); | |
1191 | Syy += cellen*pos.Y()*pos.Y(); | |
1192 | Sxy += cellen*pos.X()*pos.Y(); | |
1193 | } | |
1194 | Xc /= Ec; | |
1195 | Yc /= Ec; | |
1196 | Sxx /= Ec; | |
1197 | Syy /= Ec; | |
1198 | Sxy /= Ec; | |
1199 | Sxx -= Xc*Xc; | |
1200 | Syy -= Yc*Yc; | |
1201 | Sxy -= Xc*Yc; | |
f0897c18 | 1202 | Sxx = TMath::Abs(Sxx); |
1203 | Syy = TMath::Abs(Syy); | |
296ea9b4 | 1204 | sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0; |
1205 | sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax)); | |
1206 | sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0; | |
1207 | sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin)); | |
d595acbb | 1208 | } |
90d5b88b | 1209 | |
6bf90832 | 1210 | //________________________________________________________________________ |
296ea9b4 | 1211 | Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) const |
6bf90832 | 1212 | { |
1213 | // Calculate number of attached cells above emin. | |
1214 | ||
6bf90832 | 1215 | AliVCaloCells *cells = fEsdCells; |
1216 | if (!cells) | |
1217 | cells = fAodCells; | |
6bf90832 | 1218 | if (!cells) |
296ea9b4 | 1219 | return 0; |
6bf90832 | 1220 | |
296ea9b4 | 1221 | Int_t n = 0; |
6bf90832 | 1222 | Int_t ncells = c->GetNCells(); |
1223 | for(Int_t j=0; j<ncells; ++j) { | |
1224 | Int_t id = TMath::Abs(c->GetCellAbsId(j)); | |
1225 | Double_t cellen = cells->GetCellAmplitude(id); | |
1226 | if (cellen>=emin) | |
1227 | ++n; | |
1228 | } | |
1229 | return n; | |
1230 | } | |
1231 | ||
296ea9b4 | 1232 | //________________________________________________________________________ |
1233 | Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const | |
1234 | { | |
1235 | // Compute isolation based on tracks. | |
1236 | ||
1237 | Double_t trkIsolation = 0; | |
1238 | Double_t rad2 = radius*radius; | |
1239 | Int_t ntrks = fSelTracks->GetEntries(); | |
1240 | for(Int_t j = 0; j<ntrks; ++j) { | |
1241 | AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j)); | |
1242 | if (!track) | |
1243 | continue; | |
1244 | Float_t eta = track->Eta(); | |
1245 | Float_t phi = track->Phi(); | |
1246 | Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi); | |
1247 | Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff; | |
1248 | if(dist>rad2) | |
1249 | continue; | |
1250 | trkIsolation += track->Pt(); | |
1251 | } | |
1252 | return trkIsolation; | |
1253 | } |