]>
Commit | Line | Data |
---|---|---|
74a4d448 | 1 | #include "TChain.h" |
2 | #include "TTree.h" | |
3 | #include "TH1F.h" | |
4 | #include "TH2F.h" | |
5 | #include "TH3D.h" | |
6 | #include "TCanvas.h" | |
7 | #include "TParticle.h" | |
8 | #include "TObjArray.h" | |
9 | ||
10 | #include "AliAnalysisTask.h" | |
11 | #include "AliAnalysisManager.h" | |
12 | ||
13 | #include "AliTHn.h" | |
14 | #include "AliStack.h" | |
15 | #include "AliMCEvent.h" | |
16 | #include "AliESDEvent.h" | |
17 | #include "AliESDInputHandler.h" | |
18 | #include "AliESDtrackCuts.h" | |
19 | #include "AliCentrality.h" | |
20 | #include "AliGenHijingEventHeader.h" | |
21 | ||
22 | #include "AliAnalysisTaskEfficiencyBFPsi.h" | |
23 | ||
24 | // --------------------------------------------------------------------- | |
25 | // | |
26 | // Task for calculating the efficiency of the Balance Function | |
27 | // for single particles and pairs | |
28 | // | |
29 | // Authors: Panos Christakoglou, Michael Weber | |
30 | // | |
31 | // --------------------------------------------------------------------- | |
32 | ||
33 | ClassImp(AliAnalysisTaskEfficiencyBFPsi) | |
34 | ||
35 | //________________________________________________________________________ | |
36 | AliAnalysisTaskEfficiencyBFPsi::AliAnalysisTaskEfficiencyBFPsi(const char *name): | |
37 | AliAnalysisTaskSE(name), fESD(0), fQAList(0), fOutputList(0), | |
38 | fHistEventStats(0), fHistCentrality(0), fHistNMult(0), | |
39 | fHistGeneratedPlus(0), fHistSurvivedPlus(0), | |
40 | fHistGeneratedMinus(0), fHistSurvivedMinus(0), | |
41 | fHistGeneratedPlusPlus(0), fHistSurvivedPlusPlus(0), | |
42 | fHistGeneratedMinusMinus(0), fHistSurvivedMinusMinus(0), | |
43 | fHistGeneratedPlusMinus(0), fHistSurvivedPlusMinus(0), | |
44 | fHistGeneratedMinusPlus(0), fHistSurvivedMinusPlus(0), | |
45 | fESDtrackCuts(0), fAnalysisMode(0), | |
46 | fCentralityEstimator("V0M"), | |
47 | fCentralityPercentileMin(0.0), fCentralityPercentileMax(5.0), | |
48 | fVxMax(3.0), fVyMax(3.0), fVzMax(10.), | |
49 | fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), | |
50 | fMaxDCAxy(3.0), fMaxDCAz(3.0), | |
51 | fMinPt(0.3), fMaxPt(1.5), fMaxEta(0.8), fEtaRangeMax(0.8), | |
52 | fPtRangeMin(0.1), fPtRangeMax(5.0), fPhiRangeMin(0.0),fPhiRangeMax(6.28) { | |
53 | // Define input and output slots here | |
54 | // Input slot #0 works with a TChain | |
55 | DefineInput(0, TChain::Class()); | |
56 | // Output slot #0 id reserved by the base class for AOD | |
57 | // Output slot #1 writes into a TH1 container | |
58 | DefineOutput(1, TList::Class()); | |
59 | DefineOutput(2, TList::Class()); | |
60 | } | |
61 | ||
62 | //________________________________________________________________________ | |
63 | void AliAnalysisTaskEfficiencyBFPsi::UserCreateOutputObjects() { | |
64 | // Create histograms | |
65 | // Called once | |
66 | ||
211b716d | 67 | // global switch disabling the reference |
68 | // (to avoid "Replacing existing TH1" if several wagons are created in train) | |
69 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
70 | TH1::AddDirectory(kFALSE); | |
71 | ||
74a4d448 | 72 | fQAList = new TList(); |
73 | fQAList->SetName("QAList"); | |
74 | fQAList->SetOwner(); | |
75 | ||
76 | fOutputList = new TList(); | |
77 | fOutputList->SetName("OutputList"); | |
78 | fOutputList->SetOwner(); | |
79 | ||
80 | //Event stats. | |
81 | TString gCutName[4] = {"Total","Offline trigger", | |
82 | "Vertex","Analyzed"}; | |
83 | fHistEventStats = new TH1F("fHistEventStats", | |
84 | "Event statistics;;N_{events}", | |
85 | 4,0.5,4.5); | |
86 | for(Int_t i = 1; i <= 4; i++) | |
87 | fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data()); | |
88 | fQAList->Add(fHistEventStats); | |
89 | ||
90 | //ESD analysis | |
91 | fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events", | |
92 | 20,0.5,20.5); | |
93 | fQAList->Add(fHistCentrality); | |
94 | ||
95 | //multiplicity (good MC tracks) | |
96 | TString histName; | |
97 | histName = "fHistNMult"; | |
98 | fHistNMult = new TH1F(histName.Data(), | |
99 | ";N_{mult.}", | |
100 | 200,0,20000); | |
101 | fQAList->Add(fHistNMult); | |
102 | ||
103 | //Setting up the AliTHn | |
104 | Int_t anaSteps = 1; // analysis steps | |
105 | Int_t iBinSingle[kVariablesSingle]; // binning for track variables | |
106 | Double_t* dBinsSingle[kVariablesSingle]; // bins for track variables | |
107 | TString axisTitleSingle[kVariablesSingle]; // axis titles for track variables | |
108 | ||
109 | // two particle histograms | |
110 | Int_t iBinPair[kVariablesPair]; // binning for track variables | |
111 | Double_t* dBinsPair[kVariablesPair]; // bins for track variables | |
112 | TString axisTitlePair[kVariablesPair]; // axis titles for track variables | |
113 | ||
114 | //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (all) | |
115 | const Int_t kNPsi2Bins = 4; | |
116 | Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5}; | |
117 | iBinSingle[0] = kNPsi2Bins; | |
118 | dBinsSingle[0] = psi2Bins; | |
119 | axisTitleSingle[0] = "#phi - #Psi_{2} (a.u.)"; | |
120 | iBinPair[0] = kNPsi2Bins; | |
121 | dBinsPair[0] = psi2Bins; | |
122 | axisTitlePair[0] = "#phi - #Psi_{2} (a.u.)"; | |
123 | ||
124 | // delta eta | |
125 | const Int_t kNDeltaEtaBins = 80; | |
126 | Double_t deltaEtaBins[kNDeltaEtaBins+1]; | |
127 | for(Int_t i = 0; i < kNDeltaEtaBins+1; i++) | |
128 | deltaEtaBins[i] = -2.0 + i * 0.05; | |
129 | iBinPair[1] = kNDeltaEtaBins; | |
130 | dBinsPair[1] = deltaEtaBins; | |
131 | axisTitlePair[1] = "#Delta #eta"; | |
132 | ||
133 | // delta phi | |
134 | const Int_t kNDeltaPhiBins = 72; | |
135 | Double_t deltaPhiBins[kNDeltaPhiBins+1]; | |
136 | for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){ | |
137 | deltaPhiBins[i] = -180.0 + i * 5.; | |
138 | } | |
139 | iBinPair[2] = kNDeltaPhiBins; | |
140 | dBinsPair[2] = deltaPhiBins; | |
141 | axisTitlePair[2] = "#Delta #phi (#circ)"; | |
142 | ||
143 | // pt(trigger-associated) | |
144 | const Int_t kNPtBins = 16; | |
145 | Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.}; | |
146 | iBinSingle[1] = kNPtBins; | |
147 | dBinsSingle[1] = ptBins; | |
148 | axisTitleSingle[1] = "p_{t}^{trig.} (GeV/c)"; | |
149 | ||
150 | iBinPair[3] = kNPtBins; | |
151 | dBinsPair[3] = ptBins; | |
152 | axisTitlePair[3] = "p_{t}^{trig.} (GeV/c)"; | |
153 | ||
154 | iBinPair[4] = kNPtBins; | |
155 | dBinsPair[4] = ptBins; | |
156 | axisTitlePair[4] = "p_{t}^{assoc.} (GeV/c)"; | |
157 | ||
158 | //=============================// | |
159 | //Generated: Single particle distributions | |
160 | fHistGeneratedPlus = new AliTHn("fHistGeneratedPlus", | |
161 | "Generated positive primaries", | |
162 | anaSteps,kVariablesSingle,iBinSingle); | |
163 | for (Int_t j = 0; j < kVariablesSingle; j++) { | |
164 | fHistGeneratedPlus->SetBinLimits(j, dBinsSingle[j]); | |
165 | fHistGeneratedPlus->SetVarTitle(j, axisTitleSingle[j]); | |
166 | } | |
167 | fOutputList->Add(fHistGeneratedPlus); | |
168 | ||
169 | fHistGeneratedMinus = new AliTHn("fHistGeneratedMinus", | |
170 | "Generated positive primaries", | |
171 | anaSteps,kVariablesSingle,iBinSingle); | |
172 | for (Int_t j = 0; j < kVariablesSingle; j++) { | |
173 | fHistGeneratedMinus->SetBinLimits(j, dBinsSingle[j]); | |
174 | fHistGeneratedMinus->SetVarTitle(j, axisTitleSingle[j]); | |
175 | } | |
176 | fOutputList->Add(fHistGeneratedMinus); | |
177 | ||
178 | //Survived: Single particle distributions | |
179 | fHistSurvivedPlus = new AliTHn("fHistSurvivedPlus", | |
180 | "Survived positive primaries", | |
181 | anaSteps,kVariablesSingle,iBinSingle); | |
182 | for (Int_t j = 0; j < kVariablesSingle; j++) { | |
183 | fHistSurvivedPlus->SetBinLimits(j, dBinsSingle[j]); | |
184 | fHistSurvivedPlus->SetVarTitle(j, axisTitleSingle[j]); | |
185 | } | |
186 | fOutputList->Add(fHistSurvivedPlus); | |
187 | ||
188 | fHistSurvivedMinus = new AliTHn("fHistSurvivedMinus", | |
189 | "Survived positive primaries", | |
190 | anaSteps,kVariablesSingle,iBinSingle); | |
191 | for (Int_t j = 0; j < kVariablesSingle; j++) { | |
192 | fHistSurvivedMinus->SetBinLimits(j, dBinsSingle[j]); | |
193 | fHistSurvivedMinus->SetVarTitle(j, axisTitleSingle[j]); | |
194 | } | |
195 | fOutputList->Add(fHistSurvivedMinus); | |
196 | ||
197 | //=============================// | |
198 | //Generated-Survived: Particle pairs +- | |
199 | fHistGeneratedPlusMinus = new AliTHn("fHistGeneratedPlusMinus", | |
200 | "Generated +- primaries", | |
201 | anaSteps,kVariablesPair,iBinPair); | |
202 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
203 | fHistGeneratedPlusMinus->SetBinLimits(j, dBinsPair[j]); | |
204 | fHistGeneratedPlusMinus->SetVarTitle(j, axisTitlePair[j]); | |
205 | } | |
206 | fOutputList->Add(fHistGeneratedPlusMinus); | |
207 | ||
208 | fHistSurvivedPlusMinus = new AliTHn("fHistSurvivedPlusMinus", | |
209 | "Survived +- primaries", | |
210 | anaSteps,kVariablesPair,iBinPair); | |
211 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
212 | fHistSurvivedPlusMinus->SetBinLimits(j, dBinsPair[j]); | |
213 | fHistSurvivedPlusMinus->SetVarTitle(j, axisTitlePair[j]); | |
214 | } | |
215 | fOutputList->Add(fHistSurvivedPlusMinus); | |
216 | ||
217 | //Generated-Survived: Particle pairs -+ | |
218 | fHistGeneratedMinusPlus = new AliTHn("fHistGeneratedMinusPlus", | |
219 | "Generated -+ primaries", | |
220 | anaSteps,kVariablesPair,iBinPair); | |
221 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
222 | fHistGeneratedMinusPlus->SetBinLimits(j, dBinsPair[j]); | |
223 | fHistGeneratedMinusPlus->SetVarTitle(j, axisTitlePair[j]); | |
224 | } | |
225 | fOutputList->Add(fHistGeneratedMinusPlus); | |
226 | ||
227 | fHistSurvivedMinusPlus = new AliTHn("fHistSurvivedMinusPlus", | |
228 | "Survived -+ primaries", | |
229 | anaSteps,kVariablesPair,iBinPair); | |
230 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
231 | fHistSurvivedMinusPlus->SetBinLimits(j, dBinsPair[j]); | |
232 | fHistSurvivedMinusPlus->SetVarTitle(j, axisTitlePair[j]); | |
233 | } | |
234 | fOutputList->Add(fHistSurvivedMinusPlus); | |
235 | ||
236 | //Generated-Survived: Particle pairs ++ | |
237 | fHistGeneratedPlusPlus = new AliTHn("fHistGeneratedPlusPlus", | |
238 | "Generated ++ primaries", | |
239 | anaSteps,kVariablesPair,iBinPair); | |
240 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
241 | fHistGeneratedPlusPlus->SetBinLimits(j, dBinsPair[j]); | |
242 | fHistGeneratedPlusPlus->SetVarTitle(j, axisTitlePair[j]); | |
243 | } | |
244 | fOutputList->Add(fHistGeneratedPlusPlus); | |
245 | ||
246 | fHistSurvivedPlusPlus = new AliTHn("fHistSurvivedPlusPlus", | |
247 | "Survived ++ primaries", | |
248 | anaSteps,kVariablesPair,iBinPair); | |
249 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
250 | fHistSurvivedPlusPlus->SetBinLimits(j, dBinsPair[j]); | |
251 | fHistSurvivedPlusPlus->SetVarTitle(j, axisTitlePair[j]); | |
252 | } | |
253 | fOutputList->Add(fHistSurvivedPlusPlus); | |
254 | ||
255 | //Generated-Survived: Particle pairs -- | |
256 | fHistGeneratedMinusMinus = new AliTHn("fHistGeneratedMinusMinus", | |
257 | "Generated -- primaries", | |
258 | anaSteps,kVariablesPair,iBinPair); | |
259 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
260 | fHistGeneratedMinusMinus->SetBinLimits(j, dBinsPair[j]); | |
261 | fHistGeneratedMinusMinus->SetVarTitle(j, axisTitlePair[j]); | |
262 | } | |
263 | fOutputList->Add(fHistGeneratedMinusMinus); | |
264 | ||
265 | fHistSurvivedMinusMinus = new AliTHn("fHistSurvivedMinusMinus", | |
266 | "Survived -- primaries", | |
267 | anaSteps,kVariablesPair,iBinPair); | |
268 | for (Int_t j = 0; j < kVariablesPair; j++) { | |
269 | fHistSurvivedMinusMinus->SetBinLimits(j, dBinsPair[j]); | |
270 | fHistSurvivedMinusMinus->SetVarTitle(j, axisTitlePair[j]); | |
271 | } | |
272 | fOutputList->Add(fHistSurvivedMinusMinus); | |
273 | //=============================// | |
274 | ||
275 | fQAList->Print(); | |
276 | fOutputList->Print(); | |
277 | ||
278 | PostData(1, fQAList); | |
279 | PostData(2, fOutputList); | |
211b716d | 280 | |
281 | TH1::AddDirectory(oldStatus); | |
282 | ||
74a4d448 | 283 | } |
284 | ||
285 | //________________________________________________________________________ | |
286 | void AliAnalysisTaskEfficiencyBFPsi::UserExec(Option_t *) { | |
287 | // Main loop | |
288 | // Called for each event | |
289 | ||
290 | // Post output data. | |
291 | //ESD analysis | |
292 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); | |
293 | if (!fESD) { | |
294 | printf("ERROR: fESD not available\n"); | |
295 | return; | |
296 | } | |
297 | ||
298 | AliMCEvent* mcEvent = MCEvent(); | |
299 | if (!mcEvent) { | |
6d6eb2df | 300 | AliError("ERROR: Could not retrieve MC event"); |
74a4d448 | 301 | return; |
302 | } | |
303 | AliStack* stack = mcEvent->Stack(); | |
304 | if (!stack) { | |
6d6eb2df | 305 | AliError("ERROR: Could not retrieve MC stack"); |
74a4d448 | 306 | return; |
307 | } | |
308 | ||
309 | // arrays for 2 particle histograms | |
310 | Int_t nMCLabelCounter = 0; | |
311 | const Int_t maxMCLabelCounter = 20000; | |
312 | ||
313 | Double_t phiMinusPsi[maxMCLabelCounter]; | |
314 | Double_t eta[maxMCLabelCounter]; | |
315 | Double_t pt[maxMCLabelCounter]; | |
316 | Double_t phi[maxMCLabelCounter]; | |
317 | Int_t level[maxMCLabelCounter]; | |
318 | Int_t charge[maxMCLabelCounter]; | |
319 | ||
320 | Double_t trackVariablesSingle[kVariablesSingle]; | |
321 | Double_t trackVariablesPair[kVariablesPair]; | |
322 | ||
323 | Double_t gReactionPlane = 0.; | |
324 | AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(mcEvent)->GenEventHeader()); | |
325 | if (headerH) { | |
326 | gReactionPlane = headerH->ReactionPlaneAngle(); | |
327 | gReactionPlane *= TMath::RadToDeg(); | |
328 | } | |
329 | ||
330 | //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fESD->GetNumberOfTracks())); | |
331 | fHistEventStats->Fill(1); //all events | |
332 | ||
333 | //Centrality stuff | |
334 | AliCentrality *centrality = fESD->GetCentrality(); | |
335 | Int_t nCentrality = 0; | |
336 | nCentrality = (Int_t)(centrality->GetCentralityPercentile(fCentralityEstimator.Data())/10.); | |
337 | ||
338 | //Printf("Centrality: %lf",centrality->GetCentralityPercentile(fCentralityEstimator.Data())); | |
339 | ||
340 | if(centrality->IsEventInCentralityClass(fCentralityPercentileMin, | |
341 | fCentralityPercentileMax, | |
342 | fCentralityEstimator.Data())) { | |
343 | fHistEventStats->Fill(2); //triggered + centrality | |
344 | fHistCentrality->Fill(nCentrality+1); | |
345 | ||
346 | //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax); | |
347 | ||
348 | if(fAnalysisMode.CompareTo("TPC") == 0 ) { | |
349 | const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC(); | |
350 | if(vertex) { | |
351 | if(vertex->GetNContributors() > 0) { | |
352 | if(vertex->GetZRes() != 0) { | |
353 | fHistEventStats->Fill(3); //events with a proper vertex | |
e690d4d0 | 354 | if(TMath::Abs(vertex->GetX()) < fVxMax) { |
355 | if(TMath::Abs(vertex->GetY()) < fVyMax) { | |
356 | if(TMath::Abs(vertex->GetZ()) < fVzMax) { | |
74a4d448 | 357 | fHistEventStats->Fill(4); //analyzed events |
358 | ||
359 | Int_t nMCParticles = mcEvent->GetNumberOfTracks(); | |
360 | TArrayI labelMCArray(nMCParticles); | |
361 | ||
362 | for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) { | |
363 | AliMCParticle *mcTrack = (AliMCParticle*) mcEvent->GetTrack(iTracks); | |
364 | if (!mcTrack) { | |
6d6eb2df | 365 | AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks)); |
74a4d448 | 366 | continue; |
367 | } | |
368 | ||
369 | //exclude particles generated out of the acceptance | |
370 | Double_t vz = mcTrack->Zv(); | |
371 | if (TMath::Abs(vz) > 50.) continue; | |
372 | ||
373 | //acceptance | |
374 | if(TMath::Abs(mcTrack->Eta()) > fEtaRangeMax) | |
375 | continue; | |
376 | if((mcTrack->Pt() > fPtRangeMax)||(mcTrack->Pt() < fPtRangeMin)) | |
377 | continue; | |
d130918a | 378 | //if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin)) |
379 | //continue; | |
74a4d448 | 380 | |
381 | TParticle* particle = mcTrack->Particle(); | |
382 | if(!particle) continue; | |
383 | if(!stack->IsPhysicalPrimary(iTracks)) continue; | |
384 | ||
385 | if(iTracks <= stack->GetNprimary()) { | |
386 | Short_t gMCCharge = mcTrack->Charge(); | |
387 | Float_t firstPhi = mcTrack->Phi()*TMath::RadToDeg(); | |
388 | Float_t firstPt = mcTrack->Pt(); | |
389 | ||
390 | // Event plane (determine psi bin) | |
391 | Double_t gPsiMinusPhi = 0.; | |
392 | Double_t gPsiMinusPhiBin = -10.; | |
393 | gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane); | |
394 | //in-plane | |
395 | if((gPsiMinusPhi <= 7.5)|| | |
396 | ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5))) | |
397 | gPsiMinusPhiBin = 0.0; | |
398 | //intermediate | |
399 | else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))|| | |
400 | ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))|| | |
401 | ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))|| | |
402 | ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5))) | |
403 | gPsiMinusPhiBin = 1.0; | |
404 | //out of plane | |
405 | else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))|| | |
406 | ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5))) | |
407 | gPsiMinusPhiBin = 2.0; | |
408 | //everything else | |
409 | else | |
410 | gPsiMinusPhiBin = 3.0; | |
411 | ||
412 | trackVariablesSingle[0] = gPsiMinusPhiBin; | |
413 | trackVariablesSingle[1] = firstPt; | |
414 | ||
415 | if(gMCCharge > 0) | |
416 | fHistGeneratedPlus->Fill(trackVariablesSingle,0,1.); | |
417 | else if(gMCCharge < 0) | |
418 | fHistGeneratedMinus->Fill(trackVariablesSingle,0,1.); | |
419 | ||
420 | labelMCArray.AddAt(iTracks,nMCLabelCounter); | |
421 | if(nMCLabelCounter >= maxMCLabelCounter){ | |
422 | AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter)); | |
423 | break; | |
424 | } | |
425 | ||
426 | //fill the arrays for 2 particle analysis | |
427 | phiMinusPsi[nMCLabelCounter] = gPsiMinusPhiBin; | |
d130918a | 428 | eta[nMCLabelCounter] = particle->Eta(); |
429 | pt[nMCLabelCounter] = particle->Pt(); | |
430 | phi[nMCLabelCounter] = particle->Phi()*TMath::RadToDeg(); | |
74a4d448 | 431 | charge[nMCLabelCounter] = gMCCharge; |
432 | // findable = generated in this case! | |
433 | ||
434 | level[nMCLabelCounter] = 1; | |
435 | nMCLabelCounter += 1; | |
436 | }//primaries | |
437 | }//loop over MC particles | |
438 | ||
439 | fHistNMult->Fill(nMCLabelCounter); | |
440 | ||
441 | //ESD track loop | |
442 | Int_t nGoodTracks = fESD->GetNumberOfTracks(); | |
443 | ||
444 | TArrayI labelArray(nGoodTracks); | |
445 | Int_t labelCounter = 0; | |
446 | for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { | |
447 | AliESDtrack* track = fESD->GetTrack(iTracks); | |
448 | //AliESDtrack* track = fESDtrackCuts->GetTPCOnlyTrack(fESD,iTracks); | |
449 | if(!track) continue; | |
450 | ||
451 | AliESDtrack *tpcOnlyTrack = new AliESDtrack(); | |
452 | ||
453 | if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) { | |
454 | delete tpcOnlyTrack; | |
455 | continue; | |
456 | } | |
457 | ||
458 | Int_t label = TMath::Abs(track->GetTPCLabel()); | |
459 | if(IsLabelUsed(labelArray,label)) continue; | |
460 | labelArray.AddAt(label,labelCounter); | |
461 | labelCounter += 1; | |
fed402ec | 462 | |
74a4d448 | 463 | Int_t mcGoods = nMCLabelCounter; |
464 | for (Int_t k = 0; k < mcGoods; k++) { | |
465 | Int_t mcLabel = labelMCArray.At(k); | |
fed402ec | 466 | |
74a4d448 | 467 | if (mcLabel != TMath::Abs(label)) continue; |
468 | if(mcLabel != label) continue; | |
469 | if(label > stack->GetNtrack()) continue; | |
470 | ||
471 | TParticle *particle = stack->Particle(label); | |
472 | if(!particle) continue; | |
473 | ||
474 | //acceptance | |
475 | if(TMath::Abs(particle->Eta()) > fEtaRangeMax) | |
476 | continue; | |
477 | if((particle->Pt() > fPtRangeMax)||(particle->Pt() < fPtRangeMin)) | |
478 | continue; | |
d130918a | 479 | //if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin)) |
480 | //continue; | |
74a4d448 | 481 | |
482 | if(!stack->IsPhysicalPrimary(label)) continue; | |
483 | ||
484 | if(label <= stack->GetNprimary()) { | |
485 | Short_t gCharge = track->Charge(); | |
486 | // track cuts + analysis kinematic cuts | |
487 | if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){ | |
488 | // survived | |
489 | level[k] = 2; | |
490 | ||
491 | Float_t firstPhi = particle->Phi()*TMath::RadToDeg(); | |
492 | Float_t firstPt = particle->Pt(); | |
493 | ||
494 | // Event plane (determine psi bin) | |
495 | Double_t gPsiMinusPhi = 0.; | |
496 | Double_t gPsiMinusPhiBin = -10.; | |
497 | gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane); | |
498 | //in-plane | |
499 | if((gPsiMinusPhi <= 7.5)|| | |
500 | ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5))) | |
501 | gPsiMinusPhiBin = 0.0; | |
502 | //intermediate | |
503 | else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))|| | |
504 | ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))|| | |
505 | ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))|| | |
506 | ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5))) | |
507 | gPsiMinusPhiBin = 1.0; | |
508 | //out of plane | |
509 | else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))|| | |
510 | ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5))) | |
511 | gPsiMinusPhiBin = 2.0; | |
512 | //everything else | |
513 | else | |
514 | gPsiMinusPhiBin = 3.0; | |
515 | ||
516 | trackVariablesSingle[0] = gPsiMinusPhiBin; | |
517 | trackVariablesSingle[1] = firstPt; | |
518 | ||
519 | if(gCharge > 0) | |
520 | fHistSurvivedPlus->Fill(trackVariablesSingle,0,1.); | |
521 | else if(gCharge < 0) | |
522 | fHistSurvivedMinus->Fill(trackVariablesSingle,0,1.); | |
523 | }//track cuts | |
524 | }//primary particles | |
525 | }//loop over the generated | |
526 | }//ESD track loop | |
527 | ||
528 | labelMCArray.Reset(); | |
529 | labelArray.Reset(); | |
530 | }//Vz cut | |
531 | }//Vy cut | |
532 | }//Vx cut | |
533 | }//Vz resolution | |
534 | }//number of contributors | |
535 | }//valid vertex | |
536 | }//TPC analysis mode | |
537 | }//centrality | |
538 | ||
539 | // Here comes the 2 particle analysis | |
540 | // loop over all good MC particles | |
541 | for (Int_t i = 0; i < nMCLabelCounter ; i++) { | |
542 | Float_t firstEta = eta[i]; | |
543 | Float_t firstPhi = phi[i]; | |
544 | Float_t firstPt = pt[i]; | |
545 | Float_t gPhisMinusPsiBin = phiMinusPsi[i]; | |
546 | for (Int_t j = 0; j < nMCLabelCounter ; j++) { | |
547 | if(i == j) continue; | |
548 | ||
549 | Float_t secondEta = eta[j]; | |
550 | Float_t secondPhi = phi[j]; | |
551 | Float_t secondPt = pt[j]; | |
552 | ||
553 | trackVariablesPair[0] = gPhisMinusPsiBin; | |
554 | trackVariablesPair[1] = firstEta - secondEta; // delta eta | |
555 | trackVariablesPair[2] = firstPhi - secondPhi; // delta phi | |
556 | if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180 | |
557 | trackVariablesPair[2] -= 360.; | |
558 | if (trackVariablesPair[2] < - 180.) | |
559 | trackVariablesPair[2] += 360.; | |
560 | trackVariablesPair[3] = firstPt; // pt trigger | |
561 | trackVariablesPair[4] = secondPt; // pt | |
562 | ||
563 | //++ pairs | |
564 | if(charge[i] > 0 && charge[j] > 0 ) { | |
565 | if(level[i] > 0 && level[j] > 0) | |
566 | fHistGeneratedPlusPlus->Fill(trackVariablesPair,0,1.); | |
567 | ||
568 | if(level[i] > 1 && level[j] > 1) | |
569 | fHistSurvivedPlusPlus->Fill(trackVariablesPair,0,1.); | |
570 | } | |
571 | ||
572 | //-- pairs | |
573 | else if(charge[i] < 0 && charge[j] < 0 ) { | |
574 | if(level[i] > 0 && level[j] > 0) | |
575 | fHistGeneratedMinusMinus->Fill(trackVariablesPair,0,1.); | |
576 | ||
577 | if(level[i] > 1 && level[j] > 1) | |
578 | fHistSurvivedMinusMinus->Fill(trackVariablesPair,0,1.); | |
579 | } | |
580 | ||
581 | //+- pairs | |
582 | else if(charge[i] > 0 && charge[j] < 0 ) { | |
583 | if(level[i] > 0 && level[j] > 0) | |
584 | fHistGeneratedPlusMinus->Fill(trackVariablesPair,0,1.); | |
585 | ||
586 | if(level[i] > 1 && level[j] > 1) | |
587 | fHistSurvivedPlusMinus->Fill(trackVariablesPair,0,1.); | |
588 | } | |
589 | ||
590 | //-+ pairs | |
591 | else if(charge[i] < 0 && charge[j] > 0 ) { | |
592 | if(level[i] > 0 && level[j] > 0) | |
593 | fHistGeneratedMinusPlus->Fill(trackVariablesPair,0,1.); | |
594 | ||
595 | if(level[i] > 1 && level[j] > 1) | |
596 | fHistSurvivedMinusPlus->Fill(trackVariablesPair,0,1.); | |
597 | } | |
598 | }//second particle loop | |
599 | }//first particle loop | |
600 | ||
601 | } | |
602 | ||
603 | //________________________________________________________________________ | |
604 | void AliAnalysisTaskEfficiencyBFPsi::Terminate(Option_t *) { | |
605 | // Draw result to the screen | |
606 | // Called once at the end of the query | |
607 | ||
608 | fOutputList = dynamic_cast<TList*> (GetOutputData(1)); | |
609 | if (!fOutputList) { | |
610 | printf("ERROR: Output list not available\n"); | |
611 | return; | |
612 | } | |
613 | } | |
614 | ||
615 | //____________________________________________________________________// | |
616 | Bool_t AliAnalysisTaskEfficiencyBFPsi::IsLabelUsed(TArrayI labelArray, Int_t label) { | |
617 | //Checks if the label is used already | |
618 | Bool_t status = kFALSE; | |
619 | for(Int_t i = 0; i < labelArray.GetSize(); i++) { | |
620 | if(labelArray.At(i) == label) | |
621 | status = kTRUE; | |
622 | } | |
623 | ||
624 | return status; | |
625 | } |