]> 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
Coverity fixes:
[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
389 AliAODHeader *header = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
390 Double_t nCentrality = 0;
1cb2a06e 391
f70e6273 392 if(fUseCentrality) {
fe5628e1 393 AliCentrality *centrality = header->GetCentralityP();
1cb2a06e 394 nCentrality =centrality->GetCentralityPercentile(fCentralityEstimator.Data());
1cb2a06e 395
fe5628e1 396
397 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
398 fCentralityPercentileMax,
399 fCentralityEstimator.Data()))
400 return;
401 else {
402 fHistEventStats->Fill(2); //triggered + centrality
403 fHistCentrality->Fill(nCentrality);
404 }
405 }
406 //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
f70e6273 407
fe5628e1 408 const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
409 if(vertex) {
410 if(vertex->GetNContributors() > 0) {
411 Double32_t fCov[6];
412 vertex->GetCovarianceMatrix(fCov);
413 if(fCov[5] != 0) {
414 fHistEventStats->Fill(3); //events with a proper vertex
415 if(TMath::Abs(vertex->GetX()) < fVxMax) { // antes Xv
416 //Printf("X Vertex: %lf", vertex->GetX());
417 //Printf("Y Vertex: %lf", vertex->GetY());
418 if(TMath::Abs(vertex->GetY()) < fVyMax) { // antes Yv
419 if(TMath::Abs(vertex->GetZ()) < fVzMax) { // antes Zv
420 //Printf("Z Vertex: %lf", vertex->GetZ());
421
422 fHistEventStats->Fill(4); //analyzed events
423 fHistVz->Fill(vertex->GetZ());
424
425 //++++++++++++++++++CONTAMINATION++++++++++++++++++//
426 Int_t nGoodAODTracks = fAOD->GetNumberOfTracks();
427 Int_t nMCParticles = mcEvent->GetNumberOfTracks();
428 TArrayI labelMCArray(nMCParticles);
429
430 for(Int_t jTracks = 0; jTracks < nGoodAODTracks; jTracks++) {
431 AliAODTrack* track = fAOD->GetTrack(jTracks);
432 if(!track) continue;
433
434 if (!track->TestFilterBit(fAODTrackCutBit))
435 continue;
436
437 //acceptance
438 if(TMath::Abs(track->Eta()) > fMaxEta)
439 continue;
440 if((track->Pt() > fMaxPt)||(track->Pt() < fMinPt))
441 continue;
442
443 Double_t phiRad = track->Phi();
444
445 Int_t label = TMath::Abs(track->GetLabel());
446 if(label > nMCParticles) continue;
447 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) mcEvent->GetTrack(label);
448 Short_t gAODmcCharge = AODmcTrack->Charge();////
449 //fHistContaminationPrimaries->Fill(track->Eta(),track->Pt(),phiDeg);
450 //if (!(AODmcTrack->IsPhysicalPrimary())) {
451 //fHistContaminationSecondaries->Fill(track->Eta(),track->Pt(),phiDeg);
452 //}
731ae973
MW		 453
454 // ==============================================================================================
455 // Partial copy from AliAnalyseLeadingTrackUE::RemoveInjectedSignals:
456 // Skip tracks that come from injected signals
457 if (fInjectedSignals)
458 {
459
460 AliAODMCParticle* mother = AODmcTrack;
461
462 // find the primary mother (if not already physical primary)
463 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
464 {
465 if (((AliAODMCParticle*)mother)->GetMother() < 0)
466 {
467 mother = 0;
468 break;
469 }
470
471 mother = (AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
472 if (!mother)
473 break;
474 }
475
476
477 if (!mother)
478 {
479 AliError(Form("WARNING: No mother found for particle %d:", AODmcTrack->GetLabel()));
480 continue;
481 }
482
483 if (mother->GetLabel() >= skipParticlesAbove)
484 {
485 //AliInfo(Form("Remove particle %d (>= %d)",mother->GetLabel(),skipParticlesAbove));
486 continue;
487 }
488 }
489 // ==============================================================================================
490
fe5628e1 491 if (AODmcTrack->IsPhysicalPrimary()) {
492 if(gAODmcCharge > 0){
493 fHistContaminationPrimariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
494 }
495 if(gAODmcCharge < 0){
496 fHistContaminationPrimariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
497 }
498 }
499 else{
500 if(gAODmcCharge > 0){
501 fHistContaminationSecondariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
502 }
503 if(gAODmcCharge < 0){
504 fHistContaminationSecondariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
505 }
506 }
507 }//loop over tracks
508 //++++++++++++++++++CONTAMINATION++++++++++++++++++//
509
510 //++++++++++++++++++EFFICIENCY+++++++++++++++++++++//
511 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
512 AliAODMCParticle *mcTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
513 if (!mcTrack) {
514 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
515 continue;
516 }
517
518 //exclude particles generated out of the acceptance
519 Double_t vz = mcTrack->Zv();
520 if (TMath::Abs(vz) > 50.) continue;
521 //acceptance
522 if(TMath::Abs(mcTrack->Eta()) > fMaxEta)
523 continue;
524 if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt))
525 continue;
526
527 if(!mcTrack->IsPhysicalPrimary()) continue;
528
529 Short_t gMCCharge = mcTrack->Charge();
530 Double_t phiRad = mcTrack->Phi();
531
532 if(gMCCharge > 0)
533 fHistGeneratedEtaPtPhiPlus->Fill(mcTrack->Eta(),
534 mcTrack->Pt(),
535 phiRad);
536 else if(gMCCharge < 0)
537 fHistGeneratedEtaPtPhiMinus->Fill(mcTrack->Eta(),
538 mcTrack->Pt(),
539 phiRad);
540
541 Bool_t labelTPC = kTRUE;
542 if(labelTPC) {
543 labelMCArray.AddAt(iTracks,nMCLabelCounter);
544 if(nMCLabelCounter >= maxMCLabelCounter){
545 AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
546 break;
1cb2a06e 547 }
fe5628e1 548 //fill the arrays for 2 particle analysis
549 eta[nMCLabelCounter] = mcTrack->Eta();
550 pt[nMCLabelCounter] = mcTrack->Pt();
551 phi[nMCLabelCounter] = mcTrack->Phi();
552 charge[nMCLabelCounter] = gMCCharge;
1cb2a06e 553
fe5628e1 554 level[nMCLabelCounter] = 1;
555 nMCLabelCounter += 1;
556 }
557 }//loop over MC particles
558
559 fHistNMult->Fill(nMCLabelCounter);
560
561 //AOD track loop
562 Int_t nGoodTracks = fAOD->GetNumberOfTracks();
563 TArrayI labelArray(nGoodTracks);
564 Int_t labelCounter = 0;
565
566 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
567 AliAODTrack *trackAOD = static_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
568 if(!trackAOD) continue;
569
f70e6273 570 //track cuts
571 if (!trackAOD->TestFilterBit(fAODTrackCutBit))
572 continue;
573
574 Int_t label = TMath::Abs(trackAOD->GetLabel());
fe5628e1 575 if(IsLabelUsed(labelArray,label)) continue;
576 labelArray.AddAt(label,labelCounter);
577 labelCounter += 1;
578
fe5628e1 579 Int_t mcGoods = nMCLabelCounter;
580 for (Int_t k = 0; k < mcGoods; k++) {
581 Int_t mcLabel = labelMCArray.At(k);
1cb2a06e 582
fe5628e1 583 if (mcLabel != TMath::Abs(label)) continue;
584 if(mcLabel != label) continue;
585 if(label > trackAOD->GetLabel()) continue;
1cb2a06e 586
fe5628e1 587 //acceptance
588 if(TMath::Abs(trackAOD->Eta()) > fMaxEta)
589 continue;
590 if((trackAOD->Pt() > fMaxPt)||(trackAOD->Pt() < fMinPt))
591 continue;
1cb2a06e 592
fe5628e1 593 Short_t gCharge = trackAOD->Charge();
eb1b5eda 594 Double_t phiRad = trackAOD->Phi();
595
596 //===========================PID (so far only for electron rejection)===============================//
597 if(fElectronRejection) {
598
599 // get the electron nsigma
600 Double_t nSigma = fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kElectron);
601 fHistNSigmaTPCvsPtbeforePID->Fill(trackAOD->Pt(),nSigma);
602
603 // check only for given momentum range
604 if( trackAOD->Pt() > fElectronRejectionMinPt && trackAOD->Pt() < fElectronRejectionMaxPt ){
605
606 //look only at electron nsigma
607 if(!fElectronOnlyRejection){
608
609 //Make the decision based on the n-sigma of electrons
610 if(TMath::Abs(nSigma) < fElectronRejectionNSigma) continue;
611 }
612 else{
613
614 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kPion));
615 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kKaon));
616 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kProton));
617
618 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
619 if(TMath::Abs(nSigma) < fElectronRejectionNSigma
620 && nSigmaPions > fElectronRejectionNSigma
621 && nSigmaKaons > fElectronRejectionNSigma
622 && nSigmaProtons > fElectronRejectionNSigma ) continue;
623 }
624 }
625
626 fHistNSigmaTPCvsPtafterPID->Fill(trackAOD->Pt(),nSigma);
627
628 }
629 //===========================end of PID (so far only for electron rejection)===============================//
fe5628e1 630
631 if(TMath::Abs(trackAOD->Eta()) < fMaxEta && trackAOD->Pt() > fMinPt&&trackAOD->Pt() < fMaxPt){
f70e6273 632 level[k] = 2;
fe5628e1 633
634 if(gCharge > 0)
635 fHistSurvivedEtaPtPhiPlus->Fill(trackAOD->Eta(),
636 trackAOD->Pt(),
637 phiRad);
638 else if(gCharge < 0)
639 fHistSurvivedEtaPtPhiMinus->Fill(trackAOD->Eta(),
640 trackAOD->Pt(),
641 phiRad);
642 }//tracks
643 }//end of mcGoods
644 }//AOD track loop
645
646 labelMCArray.Reset();
647 labelArray.Reset();
648
649 }//Vz cut
650 }//Vy cut
651 }//Vx cut
652 }//Vz resolution
653 }//number of contributors
654 }//valid vertex
1cb2a06e 655
656 // Here comes the 2 particle analysis
657 // loop over all good MC particles
658 for (Int_t i = 0; i < nMCLabelCounter ; i++) {
1cb2a06e 659 // control 1D histograms (charge might be different?)
660 if(charge[i] > 0){
661 if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
f70e6273 662 if(level[i] > 1) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
1cb2a06e 663 }
664 else if(charge[i] < 0){
665 if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
f70e6273 666 if(level[i] > 1) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
1cb2a06e 667 }
668
669
670 for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
671
672 if(charge[i] > 0 && charge[j] > 0 ){
673 if(level[i] > 0 && level[j] > 0) {
674 fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 675 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 676 fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
677 }
f70e6273 678 if(level[i] > 1 && level[j] > 1) {
1cb2a06e 679 fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 680 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 681 fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
682 }
683 }
684
685 else if(charge[i] < 0 && charge[j] < 0 ){
686 if(level[i] > 0 && level[j] > 0) {
687 fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 688 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 689 fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
690 }
f70e6273 691 if(level[i] > 2 && level[j] > 1) {
1cb2a06e 692 fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 693 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 694 fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
695 }
696 }
697
698 else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
699 if(level[i] > 0 && level[j] > 0) {
700 fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 701 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 702 fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
703 }
f70e6273 704 if(level[i] > 2 && level[j] > 1) {
1cb2a06e 705 fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
f70e6273 706 if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
1cb2a06e 707 fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
708 }
709 }
710 }
711 }
1cb2a06e 712}
713
714//________________________________________________________________________
715void AliAnalysisTaskEffContBF::Terminate(Option_t *) {
716 // Draw result to the screen
717 // Called once at the end of the query
718
719 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
720 if (!fOutputList) {
721 printf("ERROR: Output list not available\n");
722 return;
723 }
724}
725
726//____________________________________________________________________//
727Bool_t AliAnalysisTaskEffContBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
728 //Checks if the label is used already
729 Bool_t status = kFALSE;
730 for(Int_t i = 0; i < labelArray.GetSize(); i++) {
731 if(labelArray.At(i) == label)
732 status = kTRUE;
733 }
734
735 return status;
736}