]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEffContBF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master_patch'
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEffContBF.cxx
CommitLineData
1cb2a06e 1#include "TChain.h"
2#include "TList.h"
3#include "TTree.h"
4#include "TH1F.h"
5#include "TH2F.h"
6#include "TH3D.h"
7#include "TCanvas.h"
8#include "TParticle.h"
9#include "TObjArray.h"
10
11#include "AliAnalysisTask.h"
12#include "AliAnalysisManager.h"
13
14#include "AliStack.h"
15#include "AliMCEvent.h"
16#include "AliAODEvent.h"
17#include "AliAODInputHandler.h"
18#include "AliAODMCParticle.h"
731ae973 19#include "AliAODMCHeader.h"
1cb2a06e 20#include "AliCentrality.h"
731ae973 21#include "AliGenEventHeader.h"
1cb2a06e 22
731ae973 23#include "AliLog.h"
1cb2a06e 24#include "AliAnalysisTaskEffContBF.h"
25
26// ---------------------------------------------------------------------
27//
28// Task for calculating the efficiency of the Balance Function
29// for single particles and pairs
30//
31// ---------------------------------------------------------------------
32
33ClassImp(AliAnalysisTaskEffContBF)
34
35//________________________________________________________________________
36AliAnalysisTaskEffContBF::AliAnalysisTaskEffContBF(const char *name)
37 : AliAnalysisTaskSE(name),
38 fAOD(0),
1cb2a06e 39 fArrayMC(0),
40 fQAList(0),
41 fOutputList(0),
42 fHistEventStats(0),
43 fHistCentrality(0),
44 fHistNMult(0),
45 fHistVz(0),
eb1b5eda 46 fHistNSigmaTPCvsPtbeforePID(0),
47 fHistNSigmaTPCvsPtafterPID(0),
27c60726 48 fHistContaminationSecondariesPlus(0),
49 fHistContaminationSecondariesMinus(0), //
50 fHistContaminationPrimariesPlus(0),
51 fHistContaminationPrimariesMinus(0), //
1cb2a06e 52 fHistGeneratedEtaPtPhiPlus(0),
53 fHistSurvivedEtaPtPhiPlus(0),
54 fHistGeneratedEtaPtPhiMinus(0),
55 fHistSurvivedEtaPtPhiMinus(0),
56 fHistGeneratedEtaPtPlusControl(0),
57 fHistSurvivedEtaPtPlusControl(0),
58 fHistGeneratedEtaPtMinusControl(0),
59 fHistSurvivedEtaPtMinusControl(0),
60 fHistGeneratedEtaPtPlusPlus(0),
61 fHistSurvivedEtaPtPlusPlus(0),
62 fHistGeneratedEtaPtMinusMinus(0),
63 fHistSurvivedEtaPtMinusMinus(0),
64 fHistGeneratedEtaPtPlusMinus(0),
65 fHistSurvivedEtaPtPlusMinus(0),
66 fHistGeneratedPhiEtaPlusPlus(0),
67 fHistSurvivedPhiEtaPlusPlus(0),
68 fHistGeneratedPhiEtaMinusMinus(0),
69 fHistSurvivedPhiEtaMinusMinus(0),
70 fHistGeneratedPhiEtaPlusMinus(0),
71 fHistSurvivedPhiEtaPlusMinus(0),
f70e6273 72 fUseCentrality(kFALSE),
1cb2a06e 73 fCentralityEstimator("V0M"),
74 fCentralityPercentileMin(0.0),
75 fCentralityPercentileMax(5.0),
eb1b5eda 76 fInjectedSignals(kFALSE),
77 fPIDResponse(0),
78 fElectronRejection(kFALSE),
79 fElectronOnlyRejection(kFALSE),
80 fElectronRejectionNSigma(-1.),
81 fElectronRejectionMinPt(0.),
82 fElectronRejectionMaxPt(1000.),
1cb2a06e 83 fVxMax(3.0),
84 fVyMax(3.0),
85 fVzMax(10.),
6d6eb2df 86 fAODTrackCutBit(128),
1cb2a06e 87 fMinNumberOfTPCClusters(80),
88 fMaxChi2PerTPCCluster(4.0),
89 fMaxDCAxy(3.0),
90 fMaxDCAz(3.0),
91 fMinPt(0.0),
92 fMaxPt(20.0),
93 fMinEta(-0.8),
94 fMaxEta(0.8),
95 fEtaRangeMin(0.0),
96 fEtaRangeMax(1.6),
97 fPtRangeMin(0.0),
98 fPtRangeMax(20.0),
1cb2a06e 99 fEtaBin(100), //=100 (BF) 16
100 fdEtaBin(64), //=64 (BF) 16
101 fPtBin(100), //=100 (BF) 36
1cb2a06e 102 fHistSurvived4EtaPtPhiPlus(0),
103 fHistSurvived8EtaPtPhiPlus(0)
104
105{
106 // Define input and output slots here
107 // Input slot #0 works with a TChain
108 DefineInput(0, TChain::Class());
109 // Output slot #0 id reserved by the base class for AOD
110 // Output slot #1 writes into a TH1 container
111 DefineOutput(1, TList::Class());
112 DefineOutput(2, TList::Class());
113}
114
115//________________________________________________________________________
116void AliAnalysisTaskEffContBF::UserCreateOutputObjects() {
117 // Create histograms
118 // Called once
119
120 fQAList = new TList();
121 fQAList->SetName("QAList");
122 fQAList->SetOwner();
123
124 fOutputList = new TList();
125 fOutputList->SetName("OutputList");
126 fOutputList->SetOwner();
127
128 //Event stats.
129 TString gCutName[4] = {"Total","Offline trigger",
130 "Vertex","Analyzed"};
131 fHistEventStats = new TH1F("fHistEventStats",
132 "Event statistics;;N_{events}",
133 4,0.5,4.5);
134 for(Int_t i = 1; i <= 4; i++)
135 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
136 fQAList->Add(fHistEventStats);
137
138 //====================================================//
139 Int_t ptBin = 36;
140 Int_t etaBin = 16;
141 Int_t phiBin = 100;
142
143 Double_t nArrayPt[37]={0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0};
144 Double_t nArrayEta[17]={-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
fe5628e1 145
146 Double_t nArrayPhi[phiBin+1];
147 for(Int_t iBin = 0; iBin <= phiBin; iBin++)
148 nArrayPhi[iBin] = iBin*TMath::TwoPi()/phiBin;
1cb2a06e 149
150 Int_t detaBin = 16;
fe5628e1 151 Int_t dphiBin = 100;
152 Double_t nArrayDPhi[dphiBin+1];
d56adb45 153 for(Int_t iBin = 0; iBin <= dphiBin; iBin++)
fe5628e1 154 nArrayDPhi[iBin] = iBin*TMath::TwoPi()/dphiBin;
1cb2a06e 155 Double_t nArrayDEta[17]={0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6};
156 //====================================================//
157
158 //AOD analysis
159 fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
fe5628e1 160 1001,-0.5,100.5);
1cb2a06e 161 fQAList->Add(fHistCentrality);
162
163 //multiplicity (good MC tracks)
164 TString histName;
165 histName = "fHistNMult";
166 fHistNMult = new TH1F(histName.Data(),
167 ";N_{mult.}",
168 200,0,20000);
169 fQAList->Add(fHistNMult);
170
171 //Vz addition+++++++++++++++++++++++++++++
172 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
173 fQAList->Add(fHistVz);
174
eb1b5eda 175 //Electron cuts -> PID QA
176 fHistNSigmaTPCvsPtbeforePID = new TH2F ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore",200, 0, 20, 200, -10, 10);
177 fQAList->Add(fHistNSigmaTPCvsPtbeforePID);
178
179 fHistNSigmaTPCvsPtafterPID = new TH2F ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter",200, 0, 20, 200, -10, 10);
180 fQAList->Add(fHistNSigmaTPCvsPtafterPID);
181
1cb2a06e 182 //Contamination for Secondaries
27c60726 183 fHistContaminationSecondariesPlus = new TH3D("fHistContaminationSecondariesPlus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
184 fOutputList->Add(fHistContaminationSecondariesPlus);
185
186 fHistContaminationSecondariesMinus = new TH3D("fHistContaminationSecondariesMinus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
187 fOutputList->Add(fHistContaminationSecondariesMinus);
1cb2a06e 188
189 //Contamination for Primaries
27c60726 190 fHistContaminationPrimariesPlus = new TH3D("fHistContaminationPrimariesPlus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
191 fOutputList->Add(fHistContaminationPrimariesPlus);
192
193 fHistContaminationPrimariesMinus = new TH3D("fHistContaminationPrimariesMinus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
194 fOutputList->Add(fHistContaminationPrimariesMinus);
1cb2a06e 195
196 //eta vs pt for MC positives
197 fHistGeneratedEtaPtPhiPlus = new TH3D("fHistGeneratedEtaPtPhiPlus",
198 "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
199 etaBin,nArrayEta, ptBin, nArrayPt,phiBin, nArrayPhi);
200 // fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,nArrayPhi);
201 fOutputList->Add(fHistGeneratedEtaPtPhiPlus);
202 fHistSurvivedEtaPtPhiPlus = new TH3D("fHistSurvivedEtaPtPhiPlus",
203 "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
204 etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
205 fOutputList->Add(fHistSurvivedEtaPtPhiPlus);
206
207 //eta vs pt for MC negatives
208 fHistGeneratedEtaPtPhiMinus = new TH3D("fHistGeneratedEtaPtPhiMinus",
209 "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
210 etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
211 fOutputList->Add(fHistGeneratedEtaPtPhiMinus);
212 fHistSurvivedEtaPtPhiMinus = new TH3D("fHistSurvivedEtaPtPhiMinus",
213 "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
214 etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
215 fOutputList->Add(fHistSurvivedEtaPtPhiMinus);
216
217 //eta vs pt for MC positives (control)
218 fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl",
219 "Generated positive primaries;#eta;p_{T} (GeV/c)",
220 etaBin,nArrayEta,ptBin,nArrayPt);
221 fOutputList->Add(fHistGeneratedEtaPtPlusControl);
222 fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl",
223 "Survived positive primaries;#eta;p_{T} (GeV/c)",
224 etaBin,nArrayEta,ptBin,nArrayPt);
225 fOutputList->Add(fHistSurvivedEtaPtPlusControl);
226
227 //eta vs pt for MC negatives (control)
228 fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl",
229 "Generated positive primaries;#eta;p_{T} (GeV/c)",
230 etaBin,nArrayEta,ptBin,nArrayPt);
231 fOutputList->Add(fHistGeneratedEtaPtMinusControl);
232 fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl",
233 "Survived positive primaries;#eta;p_{T} (GeV/c)",
234 etaBin,nArrayEta,ptBin,nArrayPt);
235 fOutputList->Add(fHistSurvivedEtaPtMinusControl);
236
237 //eta vs pt for MC ++
238 fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus",
239 "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)",
240 detaBin,nArrayDEta,ptBin,nArrayPt);
241 fOutputList->Add(fHistGeneratedEtaPtPlusPlus);
242 fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus",
243 "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)",
244 detaBin,nArrayDEta,ptBin,nArrayPt);
245 fOutputList->Add(fHistSurvivedEtaPtPlusPlus);
246
247 //eta vs pt for MC --
248 fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus",
249 "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)",
250 detaBin,nArrayDEta,ptBin,nArrayPt);
251 fOutputList->Add(fHistGeneratedEtaPtMinusMinus);
252 fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus",
253 "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)",
254 detaBin,nArrayDEta,ptBin,nArrayPt);
255 fOutputList->Add(fHistSurvivedEtaPtMinusMinus);
256
257 //eta vs pt for MC +-
258 fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus",
259 "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)",
260 detaBin,nArrayDEta,ptBin,nArrayPt);
261 fOutputList->Add(fHistGeneratedEtaPtPlusMinus);
262 fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus",
263 "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)",
264 detaBin,nArrayDEta,ptBin,nArrayPt);
265 fOutputList->Add(fHistSurvivedEtaPtPlusMinus);
266
267 //=============================//
268 //phi vs eta for MC ++
269 fHistGeneratedPhiEtaPlusPlus = new TH2F("fHistGeneratedPhiEtaPlusPlus",
270 "Generated ++ primaries;#Delta#phi",
271 dphiBin,nArrayDPhi,detaBin,nArrayDEta);
272 fOutputList->Add(fHistGeneratedPhiEtaPlusPlus);
273 fHistSurvivedPhiEtaPlusPlus = new TH2F("fHistSurvivedPhiEtaPlusPlus",
274 "Survived ++ primaries;#Delta#phi;#Delta#eta",
275 dphiBin,nArrayDPhi,detaBin,nArrayDEta);
276 fOutputList->Add(fHistSurvivedPhiEtaPlusPlus);
277
278 //phi vs eta for MC --
279 fHistGeneratedPhiEtaMinusMinus = new TH2F("fHistGeneratedPhiEtaMinusMinus",
280 "Generated -- primaries;#Delta#phi;#Delta#eta",
281 dphiBin,nArrayDPhi,detaBin,nArrayDEta);
282 fOutputList->Add(fHistGeneratedPhiEtaMinusMinus);
283 fHistSurvivedPhiEtaMinusMinus = new TH2F("fHistSurvivedPhiEtaMinusMinus",
284 "Survived -- primaries;#Delta#phi;#Delta#eta",
285 dphiBin,nArrayDPhi,detaBin,nArrayDEta);
286 fOutputList->Add(fHistSurvivedPhiEtaMinusMinus);
287
288 //phi vs eta for MC +-
289 fHistGeneratedPhiEtaPlusMinus = new TH2F("fHistGeneratedPhiEtaPlusMinus",
290 "Generated +- primaries;#Delta#phi;#Delta#eta",
291 dphiBin,nArrayDPhi,detaBin,nArrayDEta);
292 fOutputList->Add(fHistGeneratedPhiEtaPlusMinus);
293 fHistSurvivedPhiEtaPlusMinus = new TH2F("fHistSurvivedPhiEtaPlusMinus",
294 "Survived +- primaries;#Delta#phi;#Delta#eta",
295 dphiBin,nArrayDPhi,detaBin,nArrayDEta);
296 fOutputList->Add(fHistSurvivedPhiEtaPlusMinus);
297
298 fHistSurvived4EtaPtPhiPlus = new TH3F("fHistSurvived4EtaPtPhiPlus",
299 "Survived4 + primaries;#eta;p_{T} (GeV/c);#phi",
300 etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
301 fOutputList->Add(fHistSurvived4EtaPtPhiPlus);
302 fHistSurvived8EtaPtPhiPlus = new TH3F("fHistSurvived8EtaPtPhiPlus",
303 "Survived8 + primaries;#eta;p_{T} (GeV/c);#phi",
304 etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
305 fOutputList->Add(fHistSurvived8EtaPtPhiPlus);
306
307 //fQAList->Print();
308 //fOutputList->Print();
309 PostData(1, fQAList);
310 PostData(2, fOutputList);
311}
312
313//________________________________________________________________________
314void AliAnalysisTaskEffContBF::UserExec(Option_t *) {
315 // Main loop
316 // Called for each event
317
318 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
319 if (!fAOD) {
320 printf("ERROR: fAOD not available\n");
321 return;
322 }
323
324 fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
325 if (!fArrayMC)
326 AliFatal("No array of MC particles found !!!"); // MW no AliFatal use return values
327
328 AliMCEvent* mcEvent = MCEvent();
329 if (!mcEvent) {
6d6eb2df 330 AliError("ERROR: Could not retrieve MC event");
1cb2a06e 331 return;
332 }
731ae973 333
eb1b5eda 334 // PID Response task active?
335 if(fElectronRejection) {
336 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
337 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
338 }
339
731ae973
MW		 340 // ==============================================================================================
341 // Copy from AliAnalysisTaskPhiCorrelations:
342 // For productions with injected signals, figure out above which label to skip particles/tracks
343 Int_t skipParticlesAbove = 0;
344 if (fInjectedSignals)
345 {
346 AliGenEventHeader* eventHeader = 0;
347 Int_t headers = 0;
348
349 // AOD only
350 AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
d0b37020 351 if (!header){
731ae973 352 AliFatal("fInjectedSignals set but no MC header found");
d0b37020 353 return;
354 }
731ae973
MW		 355
356 headers = header->GetNCocktailHeaders();
357 eventHeader = header->GetCocktailHeader(0);
358
359
360 if (!eventHeader)
361 {
362 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
363 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
364 AliError("First event header not found. Skipping this event.");
365 return;
366 }
367
368 skipParticlesAbove = eventHeader->NProduced();
369 AliInfo(Form("Injected signals in this event (%d headers). Keeping particles/tracks of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
370 }
371 // ==============================================================================================
372
1cb2a06e 373
374 // arrays for 2 particle histograms
375 Int_t nMCLabelCounter = 0;
376 const Int_t maxMCLabelCounter = 20000;
377
378 Double_t eta[maxMCLabelCounter];
379 Double_t pt[maxMCLabelCounter];
380 Double_t phi[maxMCLabelCounter];
381 Int_t level[maxMCLabelCounter];
382 Int_t charge[maxMCLabelCounter];
383
384 //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fAOD->GetNumberOfTracks()));
385
386 fHistEventStats->Fill(1); //all events
387
388 //Centrality stuff
1cb2a06e 389 Double_t nCentrality = 0;
f70e6273 390 if(fUseCentrality) {
d3c17649 391
392 AliAODHeader *headerAOD = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
393 if (!headerAOD){
394 AliFatal("AOD header found");
395 return;
396 }
397
398 AliCentrality *centrality = headerAOD->GetCentralityP();
1cb2a06e 399 nCentrality =centrality->GetCentralityPercentile(fCentralityEstimator.Data());
1cb2a06e 400
fe5628e1 401
402 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
403 fCentralityPercentileMax,
404 fCentralityEstimator.Data()))
405 return;
406 else {
407 fHistEventStats->Fill(2); //triggered + centrality
408 fHistCentrality->Fill(nCentrality);
409 }
410 }
411 //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
f70e6273 412
fe5628e1 413 const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
414 if(vertex) {
415 if(vertex->GetNContributors() > 0) {
416 Double32_t fCov[6];
417 vertex->GetCovarianceMatrix(fCov);
418 if(fCov[5] != 0) {
419 fHistEventStats->Fill(3); //events with a proper vertex
420 if(TMath::Abs(vertex->GetX()) < fVxMax) { // antes Xv
421 //Printf("X Vertex: %lf", vertex->GetX());
422 //Printf("Y Vertex: %lf", vertex->GetY());
423 if(TMath::Abs(vertex->GetY()) < fVyMax) { // antes Yv
424 if(TMath::Abs(vertex->GetZ()) < fVzMax) { // antes Zv
425 //Printf("Z Vertex: %lf", vertex->GetZ());
426
427 fHistEventStats->Fill(4); //analyzed events
428 fHistVz->Fill(vertex->GetZ());
429
430 //++++++++++++++++++CONTAMINATION++++++++++++++++++//
431 Int_t nGoodAODTracks = fAOD->GetNumberOfTracks();
432 Int_t nMCParticles = mcEvent->GetNumberOfTracks();
433 TArrayI labelMCArray(nMCParticles);
434
435 for(Int_t jTracks = 0; jTracks < nGoodAODTracks; jTracks++) {
f15c1f69 436 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(jTracks));
437 if(!track) AliFatal("Not a standard AOD");
fe5628e1 438 if(!track) continue;
439
440 if (!track->TestFilterBit(fAODTrackCutBit))
441 continue;
442
443 //acceptance
444 if(TMath::Abs(track->Eta()) > fMaxEta)
445 continue;
446 if((track->Pt() > fMaxPt)||(track->Pt() < fMinPt))
447 continue;
448
449 Double_t phiRad = track->Phi();
450
451 Int_t label = TMath::Abs(track->GetLabel());
452 if(label > nMCParticles) continue;
453 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) mcEvent->GetTrack(label);
454 Short_t gAODmcCharge = AODmcTrack->Charge();////
455 //fHistContaminationPrimaries->Fill(track->Eta(),track->Pt(),phiDeg);
456 //if (!(AODmcTrack->IsPhysicalPrimary())) {
457 //fHistContaminationSecondaries->Fill(track->Eta(),track->Pt(),phiDeg);
458 //}
731ae973
MW		 459
460 // ==============================================================================================
461 // Partial copy from AliAnalyseLeadingTrackUE::RemoveInjectedSignals:
462 // Skip tracks that come from injected signals
463 if (fInjectedSignals)
464 {
465
466 AliAODMCParticle* mother = AODmcTrack;
467
468 // find the primary mother (if not already physical primary)
469 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
470 {
471 if (((AliAODMCParticle*)mother)->GetMother() < 0)
472 {
473 mother = 0;
474 break;
475 }
476
477 mother = (AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
478 if (!mother)
479 break;
480 }
481
482
483 if (!mother)
484 {
485 AliError(Form("WARNING: No mother found for particle %d:", AODmcTrack->GetLabel()));
486 continue;
487 }
488
489 if (mother->GetLabel() >= skipParticlesAbove)
490 {
491 //AliInfo(Form("Remove particle %d (>= %d)",mother->GetLabel(),skipParticlesAbove));
492 continue;
493 }
494 }
495 // ==============================================================================================
496
fe5628e1 497 if (AODmcTrack->IsPhysicalPrimary()) {
498 if(gAODmcCharge > 0){
499 fHistContaminationPrimariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
500 }
501 if(gAODmcCharge < 0){
502 fHistContaminationPrimariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
503 }
504 }
505 else{
506 if(gAODmcCharge > 0){
507 fHistContaminationSecondariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
508 }
509 if(gAODmcCharge < 0){
510 fHistContaminationSecondariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
511 }
512 }
513 }//loop over tracks
514 //++++++++++++++++++CONTAMINATION++++++++++++++++++//
515
516 //++++++++++++++++++EFFICIENCY+++++++++++++++++++++//
517 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
518 AliAODMCParticle *mcTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
519 if (!mcTrack) {
520 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
521 continue;
522 }
523
524 //exclude particles generated out of the acceptance
525 Double_t vz = mcTrack->Zv();
526 if (TMath::Abs(vz) > 50.) continue;
527 //acceptance
528 if(TMath::Abs(mcTrack->Eta()) > fMaxEta)
529 continue;
530 if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt))
531 continue;
532
533 if(!mcTrack->IsPhysicalPrimary()) continue;
534
535 Short_t gMCCharge = mcTrack->Charge();
536 Double_t phiRad = mcTrack->Phi();
537
538 if(gMCCharge > 0)
539 fHistGeneratedEtaPtPhiPlus->Fill(mcTrack->Eta(),
540 mcTrack->Pt(),
541 phiRad);
542 else if(gMCCharge < 0)
543 fHistGeneratedEtaPtPhiMinus->Fill(mcTrack->Eta(),
544 mcTrack->Pt(),
545 phiRad);
546
547 Bool_t labelTPC = kTRUE;
548 if(labelTPC) {
549 labelMCArray.AddAt(iTracks,nMCLabelCounter);
550 if(nMCLabelCounter >= maxMCLabelCounter){
551 AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
552 break;
1cb2a06e 553 }
fe5628e1 554 //fill the arrays for 2 particle analysis
555 eta[nMCLabelCounter] = mcTrack->Eta();
556 pt[nMCLabelCounter] = mcTrack->Pt();
557 phi[nMCLabelCounter] = mcTrack->Phi();
558 charge[nMCLabelCounter] = gMCCharge;
1cb2a06e 559
fe5628e1 560 level[nMCLabelCounter] = 1;
561 nMCLabelCounter += 1;
562 }
563 }//loop over MC particles
564
565 fHistNMult->Fill(nMCLabelCounter);
566
567 //AOD track loop
568 Int_t nGoodTracks = fAOD->GetNumberOfTracks();
569 TArrayI labelArray(nGoodTracks);
570 Int_t labelCounter = 0;
571
572 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
f15c1f69 573 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
fe5628e1 574 if(!trackAOD) continue;
575
f70e6273 576 //track cuts
577 if (!trackAOD->TestFilterBit(fAODTrackCutBit))
578 continue;
579
580 Int_t label = TMath::Abs(trackAOD->GetLabel());
fe5628e1 581 if(IsLabelUsed(labelArray,label)) continue;
582 labelArray.AddAt(label,labelCounter);
583 labelCounter += 1;
584
fe5628e1 585 Int_t mcGoods = nMCLabelCounter;
586 for (Int_t k = 0; k < mcGoods; k++) {
587 Int_t mcLabel = labelMCArray.At(k);
1cb2a06e 588
fe5628e1 589 if (mcLabel != TMath::Abs(label)) continue;
590 if(mcLabel != label) continue;
591 if(label > trackAOD->GetLabel()) continue;
1cb2a06e 592
fe5628e1 593 //acceptance
594 if(TMath::Abs(trackAOD->Eta()) > fMaxEta)
595 continue;
596 if((trackAOD->Pt() > fMaxPt)||(trackAOD->Pt() < fMinPt))
597 continue;
1cb2a06e 598
fe5628e1 599 Short_t gCharge = trackAOD->Charge();
eb1b5eda 600 Double_t phiRad = trackAOD->Phi();
601
602 //===========================PID (so far only for electron rejection)===============================//
603 if(fElectronRejection) {
604
605 // get the electron nsigma
606 Double_t nSigma = fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kElectron);
607 fHistNSigmaTPCvsPtbeforePID->Fill(trackAOD->Pt(),nSigma);
608
609 // check only for given momentum range
610 if( trackAOD->Pt() > fElectronRejectionMinPt && trackAOD->Pt() < fElectronRejectionMaxPt ){
611
612 //look only at electron nsigma
613 if(!fElectronOnlyRejection){
614
615 //Make the decision based on the n-sigma of electrons
616 if(TMath::Abs(nSigma) < fElectronRejectionNSigma) continue;
617 }
618 else{
619
620 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kPion));
621 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kKaon));
622 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kProton));
623
624 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
625 if(TMath::Abs(nSigma) < fElectronRejectionNSigma
626 && nSigmaPions > fElectronRejectionNSigma
627 && nSigmaKaons > fElectronRejectionNSigma
628 && nSigmaProtons > fElectronRejectionNSigma ) continue;
629 }
630 }
631
632 fHistNSigmaTPCvsPtafterPID->Fill(trackAOD->Pt(),nSigma);
633
634 }
635 //===========================end of PID (so far only for electron rejection)===============================//
fe5628e1 636
637 if(TMath::Abs(trackAOD->Eta()) < fMaxEta && trackAOD->Pt() > fMinPt&&trackAOD->Pt() < fMaxPt){
f70e6273 638 level[k] = 2;
fe5628e1 639
640 if(gCharge > 0)
641 fHistSurvivedEtaPtPhiPlus->Fill(trackAOD->Eta(),
642 trackAOD->Pt(),
643 phiRad);
644 else if(gCharge < 0)
645 fHistSurvivedEtaPtPhiMinus->Fill(trackAOD->Eta(),
646 trackAOD->Pt(),
647 phiRad);
648 }//tracks
649 }//end of mcGoods
650 }//AOD track loop
651
652 labelMCArray.Reset();
653 labelArray.Reset();
654
655 }//Vz cut
656 }//Vy cut
657 }//Vx cut
658 }//Vz resolution
659 }//number of contributors
660 }//valid vertex
1cb2a06e 661
662 // Here comes the 2 particle analysis
663 // loop over all good MC particles
664 for (Int_t i = 0; i < nMCLabelCounter ; i++) {
1cb2a06e 665 // control 1D histograms (charge might be different?)
666 if(charge[i] > 0){
667 if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
f70e6273 668 if(level[i] > 1) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
1cb2a06e 669 }
670 else if(charge[i] < 0){
671 if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
f70e6273 672 if(level[i] > 1) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
1cb2a06e 673 }
674
675
676 for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
677
678 if(charge[i] > 0 && charge[j] > 0 ){
679 if(level[i] > 0 && level[j] > 0) {
680 fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 681 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 682 fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
683 }
f70e6273 684 if(level[i] > 1 && level[j] > 1) {
1cb2a06e 685 fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 686 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 687 fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
688 }
689 }
690
691 else if(charge[i] < 0 && charge[j] < 0 ){
692 if(level[i] > 0 && level[j] > 0) {
693 fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 694 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 695 fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
696 }
f70e6273 697 if(level[i] > 2 && level[j] > 1) {
1cb2a06e 698 fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 699 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 700 fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
701 }
702 }
703
704 else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
705 if(level[i] > 0 && level[j] > 0) {
706 fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 707 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 708 fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
709 }
f70e6273 710 if(level[i] > 2 && level[j] > 1) {
1cb2a06e 711 fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 712 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 713 fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
714 }
715 }
716 }
717 }
1cb2a06e 718}
719
720//________________________________________________________________________
721void AliAnalysisTaskEffContBF::Terminate(Option_t *) {
722 // Draw result to the screen
723 // Called once at the end of the query
724
725 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
726 if (!fOutputList) {
727 printf("ERROR: Output list not available\n");
728 return;
729 }
730}
731
732//____________________________________________________________________//
733Bool_t AliAnalysisTaskEffContBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
734 //Checks if the label is used already
735 Bool_t status = kFALSE;
736 for(Int_t i = 0; i < labelArray.GetSize(); i++) {
737 if(labelArray.At(i) == label)
738 status = kTRUE;
739 }
740
741 return status;
742}