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