]>
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> |
717fe7de | 8 | #include <TH1F.h> |
9 | #include <TH2F.h> | |
10 | #include <TList.h> | |
11 | #include <TLorentzVector.h> | |
f5d4ab70 | 12 | #include <TNtuple.h> |
717fe7de | 13 | #include "AliAODEvent.h" |
14 | #include "AliAODVertex.h" | |
15 | #include "AliAnalysisManager.h" | |
16 | #include "AliAnalysisTaskEMCALClusterizeFast.h" | |
ea3fd2d5 | 17 | #include "AliCentrality.h" |
ea3fd2d5 | 18 | #include "AliEMCALGeoUtils.h" |
717fe7de | 19 | #include "AliESDEvent.h" |
20 | #include "AliESDVertex.h" | |
43cfaa06 | 21 | #include "AliLog.h" |
ea3fd2d5 | 22 | |
23 | ClassImp(AliAnalysisTaskEMCALPi0PbPb) | |
717fe7de | 24 | |
ea3fd2d5 | 25 | //________________________________________________________________________ |
26 | AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name) | |
27 | : AliAnalysisTaskSE(name), | |
717fe7de | 28 | fCentVar("V0M"), |
29 | fCentFrom(0), | |
30 | fCentTo(100), | |
76332037 | 31 | fVtxZMin(-10), |
32 | fVtxZMax(+10), | |
43cfaa06 | 33 | fUseQualFlag(1), |
717fe7de | 34 | fClusName(), |
f5d4ab70 | 35 | fDoNtuple(0), |
a49742b5 | 36 | fDoAfterburner(0), |
f224d35b | 37 | fAsymMax(1), |
a49742b5 | 38 | fNminCells(1), |
3c661da5 | 39 | fMinE(0.100), |
f224d35b | 40 | fMinErat(0), |
41 | fMinEcc(0), | |
d9f26424 | 42 | fNEvs(0), |
d595acbb | 43 | fGeom(0), |
717fe7de | 44 | fOutput(0), |
45 | fEsdEv(0), | |
46 | fAodEv(0), | |
43cfaa06 | 47 | fRecPoints(0), |
717fe7de | 48 | fEsdClusters(0), |
49 | fEsdCells(0), | |
50 | fAodClusters(0), | |
286b47a5 | 51 | fAodCells(0), |
fa443410 | 52 | fPtRanges(0), |
f5d4ab70 | 53 | fNtuple(0), |
fa443410 | 54 | fHCuts(0x0), |
55 | fHVertexZ(0x0), | |
76332037 | 56 | fHVertexZ2(0x0), |
fa443410 | 57 | fHCent(0x0), |
58 | fHCentQual(0x0), | |
d595acbb | 59 | fHColuRow(0x0), |
60 | fHColuRowE(0x0), | |
61 | fHCellMult(0x0), | |
62 | fHCellE(0x0), | |
63 | fHCellH(0x0), | |
6eb6260e | 64 | fHCellM(0x0), |
a49742b5 | 65 | fHCellM2(0x0), |
6eb6260e | 66 | fHClustEccentricity(0), |
fa443410 | 67 | fHClustEtaPhi(0x0), |
68 | fHClustEnergyPt(0x0), | |
69 | fHClustEnergySigma(0x0), | |
d595acbb | 70 | fHClustSigmaSigma(0x0), |
6eb6260e | 71 | fHClustNCellEnergyRatio(0x0), |
fa443410 | 72 | fHPionEtaPhi(0x0), |
73 | fHPionMggPt(0x0), | |
6eb6260e | 74 | fHPionMggAsym(0x0), |
75 | fHPionMggDgg(0x0) | |
ea3fd2d5 | 76 | { |
717fe7de | 77 | // Constructor. |
ea3fd2d5 | 78 | |
d595acbb | 79 | if (!name) |
80 | return; | |
81 | SetName(name); | |
ea3fd2d5 | 82 | DefineInput(0, TChain::Class()); |
ea3fd2d5 | 83 | DefineOutput(1, TList::Class()); |
fa443410 | 84 | fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells. " |
29b496d5 | 85 | "AOD:header,vertices,emcalCells"; |
ea3fd2d5 | 86 | } |
ea3fd2d5 | 87 | |
717fe7de | 88 | //________________________________________________________________________ |
89 | AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb() | |
90 | { | |
91 | // Destructor. | |
ea3fd2d5 | 92 | |
717fe7de | 93 | delete fOutput; fOutput = 0; |
fa443410 | 94 | delete fPtRanges; fPtRanges = 0; |
d595acbb | 95 | delete fGeom; fGeom = 0; |
96 | delete [] fHColuRow; | |
97 | delete [] fHColuRowE; | |
98 | delete [] fHCellMult; | |
ea3fd2d5 | 99 | } |
717fe7de | 100 | |
ea3fd2d5 | 101 | //________________________________________________________________________ |
102 | void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects() | |
103 | { | |
717fe7de | 104 | // Create user objects here. |
ea3fd2d5 | 105 | |
717fe7de | 106 | fOutput = new TList(); |
107 | fOutput->SetOwner(); | |
ea3fd2d5 | 108 | |
f5d4ab70 | 109 | if (fDoNtuple) { |
110 | TFile *f = OpenFile(1); | |
111 | if (f) | |
6eb6260e | 112 | fNtuple = new TNtuple(Form("nt%.0fto%.0f",fCentFrom,fCentTo),"nt", |
a49742b5 | 113 | "run:evt:cent:pt:e:emax:n:db:disp:mn:ms:chi:cpv:ecc:sig:eta:phi"); |
f5d4ab70 | 114 | } |
115 | ||
d595acbb | 116 | // histograms |
a49742b5 | 117 | TH1::SetDefaultSumw2(kTRUE); |
118 | TH2::SetDefaultSumw2(kTRUE); | |
fa443410 | 119 | fHCuts = new TH1F("hEventCuts","",4,0.5,4.5); |
120 | fHCuts->GetXaxis()->SetBinLabel(1,"All (PS)"); | |
121 | fHCuts->GetXaxis()->SetBinLabel(2,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo)); | |
122 | fHCuts->GetXaxis()->SetBinLabel(3,"QFlag"); | |
123 | fHCuts->GetXaxis()->SetBinLabel(4,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax)); | |
124 | fOutput->Add(fHCuts); | |
125 | fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25); | |
126 | fHVertexZ->SetXTitle("z [cm]"); | |
127 | fOutput->Add(fHVertexZ); | |
76332037 | 128 | fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25); |
129 | fHVertexZ2->SetXTitle("z [cm]"); | |
130 | fOutput->Add(fHVertexZ2); | |
fa443410 | 131 | fHCent = new TH1F("hCentBeforeCut","",101,-1,100); |
132 | fHCent->SetXTitle(fCentVar.Data()); | |
76332037 | 133 | fOutput->Add(fHCent); |
fa443410 | 134 | fHCentQual = new TH1F("hCentAfterCut","",101,-1,100); |
135 | fHCentQual->SetXTitle(fCentVar.Data()); | |
136 | fOutput->Add(fHCentQual); | |
90d5b88b | 137 | |
d595acbb | 138 | // histograms for cells |
139 | fHColuRow = new TH2F*[4]; | |
140 | fHColuRowE = new TH2F*[4]; | |
141 | fHCellMult = new TH1F*[4]; | |
142 | for (Int_t i = 0; i<4; ++i) { | |
143 | fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",49,0,49,25,0,25); | |
144 | fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i)); | |
145 | fHColuRow[i]->SetXTitle("col (i#eta)"); | |
146 | fHColuRow[i]->SetYTitle("row (i#phi)"); | |
147 | fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",49,0,49,25,0,25); | |
90d5b88b | 148 | fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i)); |
d595acbb | 149 | fHColuRowE[i]->SetXTitle("col (i#eta)"); |
150 | fHColuRowE[i]->SetYTitle("row (i#phi)"); | |
151 | fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000); | |
152 | fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i)); | |
90d5b88b | 153 | fHCellMult[i]->SetXTitle("# of cells"); |
d595acbb | 154 | fOutput->Add(fHColuRow[i]); |
155 | fOutput->Add(fHColuRowE[i]); | |
156 | fOutput->Add(fHCellMult[i]); | |
157 | } | |
a49742b5 | 158 | fHCellE = new TH1F("hCellE","",150,0.,15.); |
d595acbb | 159 | fHCellE->SetXTitle("E_{cell} [GeV]"); |
160 | fOutput->Add(fHCellE); | |
a49742b5 | 161 | fHCellH = new TH1F ("fHCellHighestE","",150,0.,15.); |
6eb6260e | 162 | fHCellH->SetXTitle("E^{max}_{cell} [GeV]"); |
d595acbb | 163 | fOutput->Add(fHCellH); |
a49742b5 | 164 | fHCellM = new TH1F ("fHCellMeanEperHitCell","",250,0.,2.5); |
6eb6260e | 165 | fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]"); |
166 | fOutput->Add(fHCellM); | |
a49742b5 | 167 | fHCellM2 = new TH1F ("fHCellMeanEperAllCells","",250,0.,1); |
f4ec884e | 168 | fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]"); |
a49742b5 | 169 | fOutput->Add(fHCellM2); |
90d5b88b | 170 | |
d595acbb | 171 | // histograms for clusters |
6eb6260e | 172 | fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1); |
173 | fHClustEccentricity->SetXTitle("#epsilon_{C}"); | |
174 | fOutput->Add(fHClustEccentricity); | |
a49742b5 | 175 | fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,500,1.2,2.2); |
fa443410 | 176 | fHClustEtaPhi->SetXTitle("#eta"); |
177 | fHClustEtaPhi->SetYTitle("#varphi"); | |
178 | fOutput->Add(fHClustEtaPhi); | |
179 | fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50); | |
180 | fHClustEnergyPt->SetXTitle("E [GeV]"); | |
6eb6260e | 181 | fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]"); |
fa443410 | 182 | fOutput->Add(fHClustEnergyPt); |
a49742b5 | 183 | fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50); |
184 | fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]"); | |
6eb6260e | 185 | fHClustEnergySigma->SetYTitle("E_{C} [GeV]"); |
d595acbb | 186 | fOutput->Add(fHClustEnergySigma); |
a49742b5 | 187 | fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50); |
6eb6260e | 188 | fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]"); |
189 | fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]"); | |
d595acbb | 190 | fOutput->Add(fHClustSigmaSigma); |
6eb6260e | 191 | fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05); |
192 | fHClustNCellEnergyRatio->SetXTitle("N_{cells}"); | |
193 | fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}"); | |
194 | fOutput->Add(fHClustNCellEnergyRatio); | |
90d5b88b | 195 | |
d595acbb | 196 | // histograms for pion candidates |
fa443410 | 197 | fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100,1.2,2.2); |
a49742b5 | 198 | fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}"); |
199 | fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}"); | |
fa443410 | 200 | fOutput->Add(fHPionEtaPhi); |
a49742b5 | 201 | fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0); |
fa443410 | 202 | fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
203 | fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]"); | |
204 | fOutput->Add(fHPionMggPt); | |
a49742b5 | 205 | fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1); |
fa443410 | 206 | fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
207 | fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]"); | |
208 | fOutput->Add(fHPionMggAsym); | |
a49742b5 | 209 | fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10); |
6eb6260e | 210 | fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
211 | fHPionMggDgg->SetYTitle("opening angle [grad]"); | |
212 | fOutput->Add(fHPionMggDgg); | |
78204d3d | 213 | const Int_t nbins = 20; |
214 | 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 | 215 | fPtRanges = new TAxis(nbins-1,xbins); |
216 | for (Int_t i = 0; i<=nbins; ++i) { | |
a49742b5 | 217 | fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2); |
fa443410 | 218 | fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]"); |
219 | if (i==0) | |
220 | fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0])); | |
221 | else if (i==nbins) | |
222 | fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50")); | |
223 | else | |
224 | fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i])); | |
fa443410 | 225 | fOutput->Add(fHPionInvMasses[i]); |
226 | } | |
d595acbb | 227 | |
717fe7de | 228 | PostData(1, fOutput); |
ea3fd2d5 | 229 | } |
230 | ||
231 | //________________________________________________________________________ | |
232 | void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *) | |
233 | { | |
717fe7de | 234 | // Called for each event. |
235 | ||
236 | if (!InputEvent()) | |
237 | return; | |
238 | ||
43cfaa06 | 239 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); |
240 | fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent()); | |
241 | if (fEsdEv) { | |
242 | am->LoadBranch("AliESDRun."); | |
243 | am->LoadBranch("AliESDHeader."); | |
244 | } else { | |
245 | fAodEv = dynamic_cast<AliAODEvent*>(InputEvent()); | |
246 | am->LoadBranch("header"); | |
247 | } | |
248 | ||
90d5b88b | 249 | if (!fGeom) { // set misalignment matrices (stored in first event) |
323834f0 | 250 | fGeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL"); //todo need to set name |
d595acbb | 251 | if (fEsdEv) { |
252 | for (Int_t i=0; i<4; ++i) | |
253 | fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i); | |
254 | } else { | |
255 | for (Int_t i=0; i<4; ++i) | |
256 | fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i); | |
257 | } | |
258 | fGeom->GetEMCGeometry(); | |
259 | } | |
260 | ||
286b47a5 | 261 | Int_t cut = 1; |
fa443410 | 262 | fHCuts->Fill(cut++); |
286b47a5 | 263 | |
717fe7de | 264 | const AliCentrality *centP = InputEvent()->GetCentrality(); |
265 | Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar); | |
fa443410 | 266 | fHCent->Fill(cent); |
717fe7de | 267 | if (cent<fCentFrom||cent>fCentTo) |
286b47a5 | 268 | return; |
43cfaa06 | 269 | |
fa443410 | 270 | fHCuts->Fill(cut++); |
286b47a5 | 271 | |
43cfaa06 | 272 | if (fUseQualFlag) { |
273 | if (centP->GetQuality()>0) | |
274 | return; | |
275 | } | |
286b47a5 | 276 | |
fa443410 | 277 | fHCentQual->Fill(cent); |
278 | fHCuts->Fill(cut++); | |
717fe7de | 279 | |
d9f26424 | 280 | // count number of accepted events |
281 | ++fNEvs; | |
282 | ||
43cfaa06 | 283 | if (fEsdEv) { |
fa443410 | 284 | am->LoadBranch("PrimaryVertex."); |
285 | am->LoadBranch("SPDVertex."); | |
286 | am->LoadBranch("TPCVertex."); | |
43cfaa06 | 287 | } else { |
288 | fAodEv = dynamic_cast<AliAODEvent*>(InputEvent()); | |
289 | am->LoadBranch("vertices"); | |
3c661da5 | 290 | if (!fAodEv) return; |
43cfaa06 | 291 | } |
292 | ||
293 | const AliVVertex *vertex = InputEvent()->GetPrimaryVertex(); | |
717fe7de | 294 | if (!vertex) |
295 | return; | |
296 | ||
fa443410 | 297 | fHVertexZ->Fill(vertex->GetZ()); |
286b47a5 | 298 | |
717fe7de | 299 | if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax) |
300 | return; | |
286b47a5 | 301 | |
fa443410 | 302 | fHCuts->Fill(cut++); |
76332037 | 303 | fHVertexZ2->Fill(vertex->GetZ()); |
717fe7de | 304 | |
43cfaa06 | 305 | fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used |
306 | fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached | |
307 | fEsdCells = 0; // will be set if ESD input used | |
308 | fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached | |
309 | // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode | |
310 | fAodCells = 0; // will be set if AOD input used | |
717fe7de | 311 | |
43cfaa06 | 312 | // deal with special output from AliAnalysisTaskEMCALClusterizeFast first |
313 | Bool_t clusattached = 0; | |
314 | Bool_t recalibrated = 0; | |
315 | if (1 && !fClusName.IsNull()) { | |
316 | AliAnalysisTaskEMCALClusterizeFast *cltask = 0; | |
317 | TObjArray *ts = am->GetTasks(); | |
318 | cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName)); | |
319 | if (cltask && cltask->GetClusters()) { | |
d595acbb | 320 | fRecPoints = const_cast<TObjArray*>(cltask->GetClusters()); |
43cfaa06 | 321 | clusattached = cltask->GetAttachClusters(); |
322 | if (cltask->GetCalibData()!=0) | |
323 | recalibrated = kTRUE; | |
324 | } | |
325 | } | |
326 | if (1 && AODEvent() && !fClusName.IsNull()) { | |
717fe7de | 327 | TList *l = AODEvent()->GetList(); |
328 | TClonesArray *clus = 0; | |
329 | if (l) { | |
330 | clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName)); | |
331 | fAodClusters = clus; | |
332 | } | |
ea3fd2d5 | 333 | } |
717fe7de | 334 | |
335 | if (fEsdEv) { // ESD input mode | |
43cfaa06 | 336 | if (1 && (!fRecPoints||clusattached)) { |
337 | if (!clusattached) | |
338 | am->LoadBranch("CaloClusters"); | |
717fe7de | 339 | TList *l = fEsdEv->GetList(); |
340 | TClonesArray *clus = 0; | |
341 | if (l) { | |
342 | clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters")); | |
343 | fEsdClusters = clus; | |
ea3fd2d5 | 344 | } |
345 | } | |
43cfaa06 | 346 | if (1) { |
347 | if (!recalibrated) | |
348 | am->LoadBranch("EMCALCells."); | |
717fe7de | 349 | fEsdCells = fEsdEv->GetEMCALCells(); |
350 | } | |
351 | } else if (fAodEv) { // AOD input mode | |
43cfaa06 | 352 | if (1 && (!fAodClusters || clusattached)) { |
353 | if (!clusattached) | |
354 | am->LoadBranch("caloClusters"); | |
717fe7de | 355 | TList *l = fAodEv->GetList(); |
356 | TClonesArray *clus = 0; | |
357 | if (l) { | |
358 | clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters")); | |
359 | fAodClusters = clus; | |
ea3fd2d5 | 360 | } |
361 | } | |
43cfaa06 | 362 | if (1) { |
363 | if (!recalibrated) | |
364 | am->LoadBranch("emcalCells"); | |
717fe7de | 365 | fAodCells = fAodEv->GetEMCALCells(); |
366 | } | |
367 | } else { | |
368 | AliFatal("Impossible to not have either pointer to ESD or AOD event"); | |
369 | } | |
ea3fd2d5 | 370 | |
43cfaa06 | 371 | if (1) { |
372 | AliDebug(2,Form("fRecPoints set: %p", fRecPoints)); | |
373 | AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters)); | |
374 | AliDebug(2,Form("fEsdCells set: %p", fEsdCells)); | |
375 | AliDebug(2,Form("fAodClusters set: %p", fAodClusters)); | |
376 | AliDebug(2,Form("fAodCells set: %p", fAodCells)); | |
377 | } | |
378 | ||
a49742b5 | 379 | if (fDoAfterburner) |
380 | ClusterAfterburner(); | |
6eb6260e | 381 | |
286b47a5 | 382 | FillCellHists(); |
383 | FillClusHists(); | |
384 | FillPionHists(); | |
323834f0 | 385 | FillOtherHists(); |
ea3fd2d5 | 386 | |
717fe7de | 387 | PostData(1, fOutput); |
ea3fd2d5 | 388 | } |
389 | ||
390 | //________________________________________________________________________ | |
391 | void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *) | |
392 | { | |
717fe7de | 393 | // Terminate called at the end of analysis. |
f5d4ab70 | 394 | |
395 | if (fNtuple) { | |
396 | TFile *f = OpenFile(1); | |
397 | if (f) | |
398 | fNtuple->Write(); | |
399 | } | |
d9f26424 | 400 | |
f224d35b | 401 | AliInfo(Form("\n%s: Accepted %lld events", GetName(), fNEvs)); |
ea3fd2d5 | 402 | } |
ea3fd2d5 | 403 | |
323834f0 | 404 | //________________________________________________________________________ |
405 | void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner() | |
406 | { | |
407 | // | |
408 | ||
409 | AliVCaloCells *cells = fEsdCells; | |
410 | if (!cells) | |
411 | cells = fAodCells; | |
412 | ||
413 | if (!cells) | |
414 | return; | |
415 | ||
416 | Int_t ncells = cells->GetNumberOfCells(); | |
417 | if (ncells<=0) | |
418 | return; | |
419 | ||
420 | Double_t cellMeanE = 0, cellSigE = 0; | |
421 | for (Int_t i = 0; i<ncells; ++i) { | |
422 | Double_t cellE = cells->GetAmplitude(i); | |
423 | cellMeanE += cellE; | |
424 | cellSigE += cellE*cellE; | |
425 | } | |
426 | cellMeanE /= ncells; | |
427 | cellSigE /= ncells; | |
428 | cellSigE -= cellMeanE*cellMeanE; | |
429 | if (cellSigE<0) | |
430 | cellSigE = 0; | |
431 | cellSigE = TMath::Sqrt(cellSigE / ncells); | |
432 | ||
433 | Double_t subE = cellMeanE - 7*cellSigE; | |
434 | if (subE<0) | |
435 | return; | |
436 | ||
437 | for (Int_t i = 0; i<ncells; ++i) { | |
438 | Short_t id=-1; | |
439 | Double_t amp=0,time=0; | |
440 | if (!cells->GetCell(i, id, amp, time)) | |
441 | continue; | |
442 | amp -= cellMeanE; | |
443 | if (amp<0.001) | |
444 | amp = 0; | |
445 | cells->SetCell(i, id, amp, time); | |
446 | } | |
447 | ||
448 | TObjArray *clusters = fEsdClusters; | |
449 | if (!clusters) | |
450 | clusters = fAodClusters; | |
451 | ||
452 | if (!clusters) | |
453 | return; | |
454 | ||
455 | Int_t nclus = clusters->GetEntries(); | |
456 | for (Int_t i = 0; i<nclus; ++i) { | |
457 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); | |
458 | if (!clus->IsEMCAL()) | |
459 | continue; | |
460 | Int_t nc = clus->GetNCells(); | |
461 | Double_t clusE = 0; | |
462 | UShort_t ids[100] = {0}; | |
463 | Double_t fra[100] = {0}; | |
464 | for (Int_t j = 0; j<nc; ++j) { | |
465 | Short_t id = TMath::Abs(clus->GetCellAbsId(j)); | |
466 | Double_t cen = cells->GetCellAmplitude(id); | |
467 | clusE += cen; | |
468 | if (cen>0) { | |
469 | ids[nc] = id; | |
470 | ++nc; | |
471 | } | |
472 | } | |
473 | if (clusE<=0) { | |
474 | clusters->RemoveAt(i); | |
475 | continue; | |
476 | } | |
477 | ||
478 | for (Int_t j = 0; j<nc; ++j) { | |
479 | Short_t id = ids[j]; | |
480 | Double_t cen = cells->GetCellAmplitude(id); | |
481 | fra[j] = cen/clusE; | |
482 | } | |
483 | clus->SetE(clusE); | |
484 | AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus); | |
485 | if (aodclus) { | |
486 | aodclus->Clear(""); | |
487 | aodclus->SetNCells(nc); | |
488 | aodclus->SetCellsAmplitudeFraction(fra); | |
489 | aodclus->SetCellsAbsId(ids); | |
490 | } | |
491 | } | |
492 | clusters->Compress(); | |
493 | } | |
494 | ||
286b47a5 | 495 | //________________________________________________________________________ |
496 | void AliAnalysisTaskEMCALPi0PbPb::FillCellHists() | |
497 | { | |
498 | // Fill histograms related to cell properties. | |
d595acbb | 499 | |
90d5b88b | 500 | AliVCaloCells *cells = fEsdCells; |
501 | if (!cells) | |
502 | cells = fAodCells; | |
503 | ||
504 | if (cells) { | |
505 | Int_t cellModCount[4] = {0,0,0,0}; | |
506 | Double_t cellMaxE = 0; | |
6eb6260e | 507 | Double_t cellMeanE = 0; |
90d5b88b | 508 | Int_t ncells = cells->GetNumberOfCells(); |
6eb6260e | 509 | if (ncells==0) |
510 | return; | |
90d5b88b | 511 | |
323834f0 | 512 | for (Int_t i = 0; i<ncells; ++i) { |
27422847 | 513 | Int_t absID = TMath::Abs(cells->GetCellNumber(i)); |
90d5b88b | 514 | Double_t cellE = cells->GetAmplitude(i); |
d595acbb | 515 | fHCellE->Fill(cellE); |
90d5b88b | 516 | if (cellE>cellMaxE) |
517 | cellMaxE = cellE; | |
6eb6260e | 518 | cellMeanE += cellE; |
d595acbb | 519 | |
90d5b88b | 520 | Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1; |
521 | Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta); | |
522 | if (!ret) { | |
523 | AliError(Form("Could not get cell index for %d", absID)); | |
524 | continue; | |
525 | } | |
526 | ++cellModCount[iSM]; | |
527 | Int_t iPhi=-1, iEta=-1; | |
528 | fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta); | |
529 | fHColuRow[iSM]->Fill(iEta,iPhi,1); | |
530 | fHColuRowE[iSM]->Fill(iEta,iPhi,cellE); | |
531 | } | |
532 | fHCellH->Fill(cellMaxE); | |
6eb6260e | 533 | cellMeanE /= ncells; |
534 | fHCellM->Fill(cellMeanE); | |
323834f0 | 535 | fHCellM2->Fill(cellMeanE*ncells/24/48/4); //hard-coded but there is a way to figure out from geometry |
536 | //Double_t avgE = totE / AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()/48/24; | |
537 | ||
90d5b88b | 538 | for (Int_t i=0; i<4; ++i) |
539 | fHCellMult[i]->Fill(cellModCount[i]); | |
540 | } | |
286b47a5 | 541 | } |
542 | ||
543 | //________________________________________________________________________ | |
544 | void AliAnalysisTaskEMCALPi0PbPb::FillClusHists() | |
545 | { | |
90d5b88b | 546 | // Fill histograms related to cluster properties. |
6eb6260e | 547 | Double_t clusterEcc = 0; |
548 | ||
fa443410 | 549 | Double_t vertex[3] = {0,0,0}; |
550 | InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); | |
551 | ||
90d5b88b | 552 | TObjArray *clusters = fEsdClusters; |
553 | if (!clusters) | |
554 | clusters = fAodClusters; | |
555 | ||
556 | if (clusters) { | |
557 | TLorentzVector clusterVec; | |
558 | Int_t nclus = clusters->GetEntries(); | |
fa443410 | 559 | for(Int_t i = 0; i<nclus; ++i) { |
90d5b88b | 560 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); |
fa443410 | 561 | if (!clus) |
562 | continue; | |
563 | if (!clus->IsEMCAL()) | |
564 | continue; | |
565 | clus->GetMomentum(clusterVec,vertex); | |
6eb6260e | 566 | Double_t maxAxis = 0, minAxis = 0; |
567 | GetSigma(clus,maxAxis,minAxis); | |
568 | if (maxAxis > 0) | |
569 | clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis)); | |
f224d35b | 570 | clus->SetChi2(clusterEcc); // store ecc in chi2 |
6eb6260e | 571 | fHClustEccentricity->Fill(clusterEcc); |
fa443410 | 572 | fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi()); |
573 | fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt()); | |
6eb6260e | 574 | fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E()); |
575 | fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis); | |
576 | fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E()); | |
a49742b5 | 577 | if (fNtuple) { |
f4ec884e | 578 | Float_t vals[17]; |
a49742b5 | 579 | vals[0] = InputEvent()->GetRunNumber(); |
881a35d9 | 580 | vals[1] = (((UInt_t)InputEvent()->GetOrbitNumber() << 12) | (UInt_t)InputEvent()->GetBunchCrossNumber()); |
d9f26424 | 581 | if (vals[1]<1) |
582 | vals[1] = fNEvs; | |
a49742b5 | 583 | vals[2] = InputEvent()->GetCentrality()->GetCentralityPercentileUnchecked(fCentVar); |
584 | vals[3] = clusterVec.Pt(); | |
585 | vals[4] = clusterVec.E(); | |
586 | vals[5] = GetMaxCellEnergy(clus); | |
587 | vals[6] = clus->GetNCells(); | |
588 | vals[7] = clus->GetDistanceToBadChannel(); | |
589 | vals[8] = clus->GetDispersion(); | |
590 | vals[9] = clus->GetM20(); | |
591 | vals[10] = clus->GetM02(); | |
592 | vals[11] = clus->Chi2(); | |
593 | vals[12] = clus->GetEmcCpvDistance(); | |
594 | vals[13] = clusterEcc; | |
f4ec884e | 595 | vals[14] = maxAxis; |
596 | vals[15] = clusterVec.Eta(); | |
597 | vals[16] = clusterVec.Phi(); | |
a49742b5 | 598 | fNtuple->Fill(vals); |
599 | } | |
fa443410 | 600 | } |
601 | } | |
286b47a5 | 602 | } |
603 | ||
604 | //________________________________________________________________________ | |
605 | void AliAnalysisTaskEMCALPi0PbPb::FillPionHists() | |
606 | { | |
607 | // Fill histograms related to pions. | |
286b47a5 | 608 | |
fa443410 | 609 | Double_t vertex[3] = {0,0,0}; |
610 | InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); | |
611 | ||
90d5b88b | 612 | TObjArray *clusters = fEsdClusters; |
613 | if (!clusters) | |
614 | clusters = fAodClusters; | |
fa443410 | 615 | |
90d5b88b | 616 | if (clusters) { |
617 | TLorentzVector clusterVec1; | |
618 | TLorentzVector clusterVec2; | |
619 | TLorentzVector pionVec; | |
620 | ||
621 | Int_t nclus = clusters->GetEntries(); | |
fa443410 | 622 | for (Int_t i = 0; i<nclus; ++i) { |
90d5b88b | 623 | AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i)); |
fa443410 | 624 | if (!clus1) |
625 | continue; | |
626 | if (!clus1->IsEMCAL()) | |
627 | continue; | |
3c661da5 | 628 | if (clus1->E()<fMinE) |
6eb6260e | 629 | continue; |
a49742b5 | 630 | if (clus1->GetNCells()<fNminCells) |
631 | continue; | |
f224d35b | 632 | if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat) |
633 | continue; | |
3c661da5 | 634 | if (clus1->Chi2()<fMinEcc) // eccentricity cut |
f224d35b | 635 | continue; |
fa443410 | 636 | clus1->GetMomentum(clusterVec1,vertex); |
637 | for (Int_t j = i+1; j<nclus; ++j) { | |
90d5b88b | 638 | AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j)); |
fa443410 | 639 | if (!clus2) |
640 | continue; | |
641 | if (!clus2->IsEMCAL()) | |
642 | continue; | |
3c661da5 | 643 | if (clus2->E()<fMinE) |
6eb6260e | 644 | continue; |
a49742b5 | 645 | if (clus2->GetNCells()<fNminCells) |
646 | continue; | |
f224d35b | 647 | if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat) |
648 | continue; | |
3c661da5 | 649 | if (clus2->Chi2()<fMinEcc) // eccentricity cut |
f224d35b | 650 | continue; |
fa443410 | 651 | clus2->GetMomentum(clusterVec2,vertex); |
652 | pionVec = clusterVec1 + clusterVec2; | |
653 | Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E(); | |
6eb6260e | 654 | Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect()); |
d595acbb | 655 | if (pionZgg < fAsymMax) { |
656 | fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi()); | |
657 | fHPionMggPt->Fill(pionVec.M(),pionVec.Pt()); | |
658 | fHPionMggAsym->Fill(pionVec.M(),pionZgg); | |
6eb6260e | 659 | fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle); |
d595acbb | 660 | Int_t bin = fPtRanges->FindBin(pionVec.Pt()); |
661 | fHPionInvMasses[bin]->Fill(pionVec.M()); | |
662 | } | |
fa443410 | 663 | } |
664 | } | |
90d5b88b | 665 | } |
fa443410 | 666 | } |
d595acbb | 667 | |
6eb6260e | 668 | //________________________________________________________________________ |
323834f0 | 669 | void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists() |
6eb6260e | 670 | { |
323834f0 | 671 | // Fill histograms related to cell properties. |
6eb6260e | 672 | } |
673 | ||
d595acbb | 674 | //________________________________________________________________________ |
90d5b88b | 675 | Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster) |
d595acbb | 676 | { |
90d5b88b | 677 | // Get maximum energy of attached cell. |
678 | ||
679 | Double_t maxe = 0; | |
90d5b88b | 680 | Int_t ncells = cluster->GetNCells(); |
681 | if (fEsdCells) { | |
682 | for (Int_t i=0; i<ncells; i++) { | |
27422847 | 683 | Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); |
90d5b88b | 684 | if (e>maxe) |
685 | maxe = e; | |
686 | } | |
687 | } else { | |
688 | for (Int_t i=0; i<ncells; i++) { | |
27422847 | 689 | Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); |
90d5b88b | 690 | if (e>maxe) |
691 | maxe = e; | |
692 | } | |
693 | } | |
6eb6260e | 694 | return maxe; |
d595acbb | 695 | } |
696 | ||
697 | //________________________________________________________________________ | |
6eb6260e | 698 | void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *cluster, Double_t& sigmaMax, Double_t &sigmaMin) |
d595acbb | 699 | { |
90d5b88b | 700 | // Calculate the (E) weighted variance along the longer (eigen) axis. |
701 | ||
6eb6260e | 702 | sigmaMax = 0; // cluster variance along its longer axis |
703 | sigmaMin = 0; // cluster variance along its shorter axis | |
90d5b88b | 704 | Double_t Ec = cluster->E(); // cluster energy |
705 | Double_t Xc = 0 ; // cluster first moment along X | |
706 | Double_t Yc = 0 ; // cluster first moment along Y | |
707 | Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2) | |
708 | Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2) | |
709 | Double_t Syy = 0 ; // cluster covariance^2 | |
710 | ||
711 | AliVCaloCells *cells = fEsdCells; | |
712 | if (!cells) | |
713 | cells = fAodCells; | |
714 | ||
6eb6260e | 715 | if (!cells) |
716 | return; | |
717 | ||
718 | Int_t ncells = cluster->GetNCells(); | |
719 | if (ncells==1) | |
720 | return; | |
721 | ||
722 | TVector3 pos; | |
723 | for(Int_t j=0; j<ncells; ++j) { | |
27422847 | 724 | Int_t id = TMath::Abs(cluster->GetCellAbsId(j)); |
725 | fGeom->GetGlobal(id,pos); | |
6eb6260e | 726 | Double_t cellen = cells->GetCellAmplitude(id); |
727 | Xc += cellen*pos.X(); | |
728 | Yc += cellen*pos.Y(); | |
729 | Sxx += cellen*pos.X()*pos.X(); | |
730 | Syy += cellen*pos.Y()*pos.Y(); | |
731 | Sxy += cellen*pos.X()*pos.Y(); | |
732 | } | |
733 | Xc /= Ec; | |
734 | Yc /= Ec; | |
735 | Sxx /= Ec; | |
736 | Syy /= Ec; | |
737 | Sxy /= Ec; | |
738 | Sxx -= Xc*Xc; | |
739 | Syy -= Yc*Yc; | |
740 | Sxy -= Xc*Yc; | |
741 | sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0; | |
742 | sigmaMax = TMath::Sqrt(sigmaMax); | |
743 | sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0; | |
744 | sigmaMin = TMath::Sqrt(sigmaMin); | |
d595acbb | 745 | } |
90d5b88b | 746 |