]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskMinijet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correct warnings
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskMinijet.cxx
CommitLineData
117f99e3 1#include <TChain.h>
2#include <TList.h>
3
4#include <TTree.h>
5#include <TH1F.h>
6#include <TH2F.h>
7099c89a 7#include <THnSparse.h>
117f99e3 8#include <TProfile.h>
9#include <TCanvas.h>
10#include "TRandom.h"
11
12#include "AliVEvent.h"
13#include "AliVParticle.h"
14
15#include "AliESDEvent.h"
16#include "AliESDtrack.h"
17#include "AliMultiplicity.h"
18#include "AliESDtrackCuts.h"
19
20#include "AliAODEvent.h"
21#include "AliAODTrack.h"
b7c0438e 22#include "AliAODMCParticle.h"
7099c89a 23#include "AliAODMCHeader.h"
117f99e3 24
25#include "AliStack.h"
26#include "AliMCEvent.h"
27#include "AliMCParticle.h"
28#include "AliGenEventHeader.h"
29
30#include "AliAnalysisManager.h"
7099c89a 31#include "AliInputEventHandler.h"
32
33#include <vector>
34#include <algorithm>
35#include <functional>
36using namespace std;
117f99e3 37
38#include "AliAnalysisTaskMinijet.h"
39
f2565a6b 40// Analysis task for two-particle correlations using all particles over pt threshold
41// pt_trig threshold for trigger particle (event axis) and pt_assoc for possible associated particles.
42// Extract mini-jet yield and fragmentation properties via Delta-Phi histograms of these correlations
43// post processing of analysis output via macro plot3and2Gaus.C
44// Can use ESD or AOD, reconstructed and Monte Carlo data as input
45// Author: Eva Sicking
46
117f99e3 47
48ClassImp(AliAnalysisTaskMinijet)
49
50//________________________________________________________________________
b9ad6f04 51 AliAnalysisTaskMinijet::AliAnalysisTaskMinijet(const char *name)
52 : AliAnalysisTaskSE(name),
53 fUseMC(kFALSE),
54 fMcOnly(kFALSE),
7099c89a 55 fBSign(0),
56 fAnalysePrimOnly(kFALSE),// not used
57 fPtMin(0.2),
58 fPtMax(50.0),
b9ad6f04 59 fCuts(0),
b9ad6f04 60 fTriggerPtCut(0.7),
61 fAssociatePtCut(0.4),
b9ad6f04 62 fMode(0),
63 fVertexZCut(10.),
64 fEtaCut(0.9),
9ed461df 65 fEtaCutSeed(0.9),
66 fSelectParticles(1),
67 fSelectParticlesAssoc(1),
b9ad6f04 68 fESDEvent(0),
69 fAODEvent(0),
9ed461df 70 fNMcPrimAccept(0),
7099c89a 71 fNRecAccept(0),
9ed461df 72 fVzEvent(0),
7099c89a 73 fMeanPtRec(0),
74 fLeadingPtRec(0),
b9ad6f04 75 fHists(0),
7099c89a 76 fStep(0),
b9ad6f04 77 fHistPt(0),
78 fHistPtMC(0),
9ed461df 79 fNmcNch(0),
f2565a6b 80 fPNmcNch(0),
b9ad6f04 81 fChargedPi0(0)
117f99e3 82{
7099c89a 83
f2565a6b 84 //Constructor
85
7099c89a 86 for(Int_t i = 0;i< 6;i++){
87 fMapSingleTrig[i] = 0;
88 fMapPair[i] = 0;
89 fMapEvent[i] = 0;
90 fMapAll[i] = 0;
91
b9ad6f04 92 fVertexZ[i] = 0;
7099c89a 93
94 fNcharge[i] = 0;
b9ad6f04 95 fPt[i] = 0;
96 fEta[i] = 0;
97 fPhi[i] = 0;
9ed461df 98 fDcaXY[i] = 0;
99 fDcaZ[i] = 0;
b9ad6f04 100
101 fPtSeed[i] = 0;
102 fEtaSeed[i] = 0;
103 fPhiSeed[i] = 0;
104
9ed461df 105 fPtOthers[i] = 0;
106 fEtaOthers[i] = 0;
107 fPhiOthers[i] = 0;
108 fPtEtaOthers[i] = 0;
b18f7dfb 109
b9ad6f04 110 fPhiEta[i] = 0;
111
112 fDPhiDEtaEventAxis[i] = 0;
9ed461df 113 fDPhiDEtaEventAxisSeeds[i]= 0;
b9ad6f04 114 fTriggerNch[i] = 0;
7099c89a 115
9ed461df 116 fTriggerNchSeeds[i] = 0;
b9ad6f04 117 fTriggerTracklet[i] = 0;
117f99e3 118
b9ad6f04 119 fNch07Nch[i] = 0;
7099c89a 120 fPNch07Nch[i] = 0;
121
b9ad6f04 122 fNch07Tracklet[i] = 0;
123 fNchTracklet[i] = 0;
7099c89a 124 fPNch07Tracklet[i] = 0;
125 fDPhiEventAxis[i] = 0;
117f99e3 126 }
127 DefineOutput(1, TList::Class());
128}
129
130//________________________________________________________________________
131AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet()
132{
133 // Destructor
134
135 if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists;
136}
137
138//________________________________________________________________________
139void AliAnalysisTaskMinijet::UserCreateOutputObjects()
140{
141 // Create histograms
142 // Called once
b9ad6f04 143 if(fDebug) Printf("In User Create Output Objects.");
117f99e3 144
7099c89a 145 fStep = new TH1F("fStep", "fStep", 10, -0.5, 9.5);
9ed461df 146 fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 150, 0.1, 3.1);
117f99e3 147 fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
148 fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
149 fHistPt->SetMarkerStyle(kFullCircle);
117f99e3 150
151 if (fUseMC) {
9ed461df 152 fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 150, 0.1, 3.1);
117f99e3 153 fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)");
154 fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
155 fHistPtMC->SetMarkerStyle(kFullCircle);
7099c89a 156
9ed461df 157 fNmcNch = new TH2F("fNmcNch", "fNmcNch", 100,-0.5,99.5,100,-0.5,99.5);
f2565a6b 158 fPNmcNch = new TProfile("pNmcNch", "pNmcNch", 100,-0.5,99.5);
9ed461df 159
117f99e3 160 }
161
b9ad6f04 162 fChargedPi0 = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5);
7099c89a 163
164
165 //----------------------
166 //bins for pt in THnSpare
167 Double_t ptMin = 0.0, ptMax = 100.;
168 Int_t nPtBins = 39;
169 Double_t binsPt[] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8,
170 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0,
171 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
172
173 // Int_t nPtBins = 100;
174 // Double_t ptMin = 1.e-2, ptMax = 100.;
175 // Double_t *binsPt = 0;
176 // binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
177
178
179 Double_t ptMin2 = 0.0, ptMax2 = 100.;
180 Int_t nPtBins2 = 9;
181 Double_t binsPt2[] = {0.1, 0.4, 0.7, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0};
182
183 // Int_t nPtBins2 = 10;
184 // Double_t ptMin2 = 0.4, ptMax2 = 100.;
185 // Double_t *binsPt2 = 0;
186 // binsPt2 = CreateLogAxis(nPtBins2,ptMin2,ptMax2);
187
188 //3 dim matrix
189 Int_t binsEffHisto[3] = {Int_t(fEtaCut*20), nPtBins, 150 };
190 Double_t minEffHisto[3] = {-fEtaCut, ptMin, -0.5 };
191 Double_t maxEffHisto[3] = {fEtaCut, ptMax, 149.5 };
192
193 //5 dim matrix
194 Int_t binsEffHisto5[6] = { nPtBins2, nPtBins2, Int_t(fEtaCut*10), 180, 150 , 2 };
195 Double_t minEffHisto5[6] = { ptMin2, ptMin2, -2*fEtaCut, -0.5*TMath::Pi(), -0.5 , -0.5 };
196 Double_t maxEffHisto5[6] = { ptMax2, ptMax2, 2*fEtaCut, 1.5*TMath::Pi(), 149.5 , 1.5 };
197
b9ad6f04 198
7099c89a 199 //4 dim matrix
200 Int_t binsEvent[4] = { 150, 20, 50, nPtBins };
201 Double_t minEvent[4] = { -0.5, -10, 0, ptMin };
202 Double_t maxEvent[4] = { 149.5, 10, 10, ptMax };
203
204 //3 dim matrix
205 Int_t binsAll[3] = {Int_t(fEtaCut*20), nPtBins2, 150 };
206 Double_t minAll[3] = {-fEtaCut, ptMin2, -0.5 };
207 Double_t maxAll[3] = {fEtaCut, ptMax2, 149.5 };
208
209 //--------------------
210 TString dataType[2] ={"ESD", "AOD"};
211 TString labels[6]={Form("%sAllAllMcNmc",dataType[fMode].Data()),
212 Form("%sTrigAllMcNmc",dataType[fMode].Data()),
213 Form("%sTrigAllMcNrec",dataType[fMode].Data()),
214 Form("%sTrigVtxMcNrec",dataType[fMode].Data()),
215 Form("%sTrigVtxRecMcPropNrec",dataType[fMode].Data()),
216 Form("%sTrigVtxRecNrec",dataType[fMode].Data())};
217
218
219 for(Int_t i=0;i<6;i++){
220
221 fMapSingleTrig[i] = new THnSparseD(Form("fMapSingleTrig%s", labels[i].Data()),"eta:pt:Nrec",
222 3,binsEffHisto,minEffHisto,maxEffHisto);
223 fMapSingleTrig[i]->SetBinEdges(1,binsPt);
224 fMapSingleTrig[i]->GetAxis(0)->SetTitle("#eta");
225 fMapSingleTrig[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
226 fMapSingleTrig[i]->GetAxis(2)->SetTitle("N_{rec}");
227 fMapSingleTrig[i]->Sumw2();
228
229 fMapPair[i] = new THnSparseD(Form("fMapPair%s", labels[i].Data()),"pt_trig:pt_assoc:DeltaEta:DeltaPhi:Nrec:likesign",
230 6,binsEffHisto5,minEffHisto5,maxEffHisto5);
231 fMapPair[i]->SetBinEdges(0,binsPt2);
232 fMapPair[i]->SetBinEdges(1,binsPt2);
233 fMapPair[i]->GetAxis(0)->SetTitle("p_{T, trig} (GeV/c)");
234 fMapPair[i]->GetAxis(1)->SetTitle("p_{T, assoc} (GeV/c)");
235 fMapPair[i]->GetAxis(2)->SetTitle("#Delta #eta");
236 fMapPair[i]->GetAxis(3)->SetTitle("#Delta #phi");
237 fMapPair[i]->GetAxis(4)->SetTitle("N_{rec}");
238 fMapPair[i]->GetAxis(5)->SetTitle("Like-sign or Unlike-sign");
239 fMapPair[i]->Sumw2();
240
241
242 fMapEvent[i] = new THnSparseD(Form("fMapEvent%s", labels[i].Data()),"Nrec:vertexZ:meanPt:leadingPt",
243 4,binsEvent,minEvent,maxEvent);
244 fMapEvent[i]->GetAxis(0)->SetTitle("N_{rec}");
245 fMapEvent[i]->GetAxis(1)->SetTitle("z_{vertex} (cm)");
246 fMapEvent[i]->GetAxis(2)->SetTitle("meanPt");
247 fMapEvent[i]->SetBinEdges(3,binsPt);
248 fMapEvent[i]->GetAxis(3)->SetTitle("leadingPt");
249 fMapEvent[i]->Sumw2();
250
251 fMapAll[i] = new THnSparseD(Form("fMapAll%s", labels[i].Data()),"eta:pt:Nrec",
252 3,binsAll,minAll,maxAll);
253 fMapAll[i]->SetBinEdges(1,binsPt2);
254 fMapAll[i]->GetAxis(0)->SetTitle("#eta");
255 fMapAll[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
256 fMapAll[i]->GetAxis(2)->SetTitle("N_{rec}");
257 fMapAll[i]->Sumw2();
117f99e3 258
117f99e3 259
260 fVertexZ[i] = new TH1F(Form("fVertexZ%s",labels[i].Data()),
b9ad6f04 261 Form("fVertexZ%s",labels[i].Data()) ,
262 200, -10., 10.);
117f99e3 263 fPt[i] = new TH1F(Form("fPt%s",labels[i].Data()),
264 Form("fPt%s",labels[i].Data()) ,
9ed461df 265 200, 0., 50);
7099c89a 266 fNcharge[i] = new TH1F(Form("fNcharge%s",labels[i].Data()),
267 Form("fNcharge%s",labels[i].Data()) ,
268 250, -0.5, 249.5);
117f99e3 269 fEta[i] = new TH1F (Form("fEta%s",labels[i].Data()),
b9ad6f04 270 Form("fEta%s",labels[i].Data()) ,
9ed461df 271 100, -2., 2);
117f99e3 272 fPhi[i] = new TH1F(Form("fPhi%s",labels[i].Data()),
b9ad6f04 273 Form("fPhi%s",labels[i].Data()) ,
274 360, 0.,2*TMath::Pi());
9ed461df 275 fDcaXY[i] = new TH1F(Form("fDcaXY%s",labels[i].Data()),
276 Form("fDcaXY%s",labels[i].Data()) ,
277 200, -3,3);
278 fDcaZ[i] = new TH1F(Form("fDcaZ%s",labels[i].Data()),
279 Form("fDcaZ%s",labels[i].Data()) ,
280 200, -10,10);
b9ad6f04 281 fPtSeed[i] = new TH1F(Form("fPSeedt%s",labels[i].Data()),
7099c89a 282 Form("fPtSeed%s",labels[i].Data()) ,
283 500, 0., 50);
b9ad6f04 284 fEtaSeed[i] = new TH1F (Form("fEtaSeed%s",labels[i].Data()),
7099c89a 285 Form("fEtaSeed%s",labels[i].Data()) ,
286 100, -2., 2);
b9ad6f04 287 fPhiSeed[i] = new TH1F(Form("fPhiSeed%s",labels[i].Data()),
7099c89a 288 Form("fPhiSeed%s",labels[i].Data()) ,
289 360, 0.,2*TMath::Pi());
b9ad6f04 290
b18f7dfb 291 fPtOthers[i] = new TH1F(Form("fPOtherst%s",labels[i].Data()),
7099c89a 292 Form("fPtOthers%s",labels[i].Data()) ,
293 500, 0., 50);
b18f7dfb 294 fEtaOthers[i] = new TH1F (Form("fEtaOthers%s",labels[i].Data()),
7099c89a 295 Form("fEtaOthers%s",labels[i].Data()) ,
296 100, -2., 2);
b18f7dfb 297 fPhiOthers[i] = new TH1F(Form("fPhiOthers%s",labels[i].Data()),
7099c89a 298 Form("fPhiOthers%s",labels[i].Data()) ,
299 360, 0.,2*TMath::Pi());
b18f7dfb 300 fPtEtaOthers[i] = new TH2F(Form("fPtEtaOthers%s",labels[i].Data()),
301 Form("fPtEtaOthers%s",labels[i].Data()) ,
302 500, 0., 50, 100, -1., 1);
b9ad6f04 303 fPhiEta[i] = new TH2F(Form("fPhiEta%s",labels[i].Data()),
304 Form("fPhiEta%s",labels[i].Data()) ,
305 180, 0., 2*TMath::Pi(), 100, -1.,1.);
117f99e3 306 fDPhiDEtaEventAxis[i] = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()),
b9ad6f04 307 Form("fDPhiDEtaEventAxis%s",labels[i].Data()) ,
7099c89a 308 180, -0.5* TMath::Pi(), 1.5*TMath::Pi(), 9, -1.8,1.8);
117f99e3 309 fTriggerNch[i] = new TH1F(Form("fTriggerNch%s",labels[i].Data()),
b9ad6f04 310 Form("fTriggerNch%s",labels[i].Data()) ,
311 250, -0.5, 249.5);
9ed461df 312 fTriggerNchSeeds[i] = new TH2F(Form("fTriggerNchSeeds%s",labels[i].Data()),
313 Form("fTriggerNchSeeds%s",labels[i].Data()) ,
314 250, -0.5, 249.5, 100, -0.5, 99.5);
117f99e3 315 fTriggerTracklet[i] = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()),
b9ad6f04 316 Form("fTriggerTracklet%s",labels[i].Data()) ,
317 250, -0.5, 249.5);
117f99e3 318 fNch07Nch[i] = new TH2F(Form("fNch07Nch%s",labels[i].Data()),
b9ad6f04 319 Form("fNch07Nch%s",labels[i].Data()) ,
320 250, -2.5, 247.5,250, -2.5, 247.5);
f2565a6b 321 fPNch07Nch[i] = new TProfile(Form("pNch07Nch%s",labels[i].Data()),
7099c89a 322 Form("pNch07Nch%s",labels[i].Data()) ,
323 250, -2.5, 247.5);
117f99e3 324 fNch07Tracklet[i] = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()),
b9ad6f04 325 Form("fNch07Tracklet%s",labels[i].Data()) ,
326 250, -2.5, 247.5,250, -2.5, 247.5);
117f99e3 327 fNchTracklet[i] = new TH2F(Form("fNchTracklet%s",labels[i].Data()),
b9ad6f04 328 Form("fNchTracklet%s",labels[i].Data()) ,
329 250, -2.5, 247.5,250, -2.5, 247.5);
f2565a6b 330 fPNch07Tracklet[i] = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()),
7099c89a 331 Form("pNch07Tracklet%s",labels[i].Data()) ,
332 250, -2.5, 247.5);
9ed461df 333 fDPhiEventAxis[i] = new TH1F(Form("fDPhiEventAxis%s",
334 labels[i].Data()),
335 Form("fDPhiEventAxis%s",
336 labels[i].Data()) ,
337 180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
338
117f99e3 339 }
340
341 fHists = new TList();
342 fHists->SetOwner();
343
7099c89a 344 fHists->Add(fStep);
117f99e3 345 fHists->Add(fHistPt);
7099c89a 346
9ed461df 347 if(fUseMC){
348 fHists->Add(fHistPtMC);
349 fHists->Add(fNmcNch);
f2565a6b 350 fHists->Add(fPNmcNch);
9ed461df 351 }
b9ad6f04 352 fHists->Add(fChargedPi0);
7099c89a 353
354 for(Int_t i=0;i<6;i++){
355 fHists->Add(fMapSingleTrig[i]);
356 fHists->Add(fMapPair[i]);
357 fHists->Add(fMapEvent[i]);
358 fHists->Add(fMapAll[i]);
117f99e3 359 fHists->Add(fVertexZ[i]);
360 fHists->Add(fPt[i]);
7099c89a 361 fHists->Add(fNcharge[i]);
117f99e3 362 fHists->Add(fEta[i]);
363 fHists->Add(fPhi[i]);
9ed461df 364 fHists->Add(fDcaXY[i]);
365 fHists->Add(fDcaZ[i]);
b9ad6f04 366 fHists->Add(fPtSeed[i]);
367 fHists->Add(fEtaSeed[i]);
368 fHists->Add(fPhiSeed[i]);
b18f7dfb 369 fHists->Add(fPtOthers[i]);
370 fHists->Add(fEtaOthers[i]);
371 fHists->Add(fPhiOthers[i]);
372 fHists->Add(fPtEtaOthers[i]);
b9ad6f04 373 fHists->Add(fPhiEta[i]);
117f99e3 374 fHists->Add(fDPhiDEtaEventAxis[i]);
375 fHists->Add(fTriggerNch[i]);
9ed461df 376 fHists->Add(fTriggerNchSeeds[i]);
117f99e3 377 fHists->Add(fTriggerTracklet[i]);
378 fHists->Add(fNch07Nch[i]);
f2565a6b 379 fHists->Add(fPNch07Nch[i]);
117f99e3 380 fHists->Add(fNch07Tracklet[i]);
381 fHists->Add(fNchTracklet[i]);
f2565a6b 382 fHists->Add(fPNch07Tracklet[i]);
9ed461df 383 fHists->Add(fDPhiEventAxis[i]);
117f99e3 384 }
7099c89a 385
117f99e3 386 PostData(1, fHists);
387
388}
389
390//________________________________________________________________________
391void AliAnalysisTaskMinijet::UserExec(Option_t *)
392{
7099c89a 393 // Main loop, called for each event
394 // Kinematics-only, ESD and AOD can be processed.
395 // Data is read (LoopESD, LoopAOD...) and then analysed (Analyse).
396 // - in case of MC with full detector simulation, all correction steps(0-5) can be processed
397 // - for Data, only step 5 is performed
398 // - for kinematics-only, only step 0 is processed
399 // step 5 = Triggered events, reconstructed accepted vertex, reconstructed tracks, reconstructed multiplicity,
400 // step 4 = Triggered events, reconstructed accepted vertex, reconstructed tracks with MC properties, reconstructed multiplicity
401 // step 3 = Triggered events, reconstructed accepted vertex, mc primary particles, reconstructed multiplicity,
402 // step 2 = Triggered events, all mc primary particles, reconstructed multiplicity
403 // step 1 = Triggered events, all mc primary particles, true multiplicity
7099c89a 404 // step 0 = All events, all mc primary particles, true multiplicity
405
406 if(fDebug) Printf("UserExec: Event starts");
9ed461df 407
117f99e3 408 AliVEvent *event = InputEvent();
409 if (!event) {
b9ad6f04 410 Error("UserExec", "Could not retrieve event");
411 return;
117f99e3 412 }
7099c89a 413 fBSign= event->GetMagneticField();
117f99e3 414
415 //get events, either ESD or AOD
416 fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
417 fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
418
7099c89a 419 vector<Float_t> pt;
420 vector<Float_t> eta;
421 vector<Float_t> phi;
422 vector<Short_t> charge;
423 vector<Float_t> px;
424 vector<Float_t> py;
425 vector<Float_t> pz;
426 vector<Float_t> theta;
427
117f99e3 428
117f99e3 429 //number of accepted tracks and tracklets
430 Int_t ntracks = 0; //return value for reading functions for ESD and AOD
7099c89a 431 //Int_t ntracksRemove = 0; //return value for reading functions for ESD and AOD
432 vector<Int_t> nTracksTracklets; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc,
433 //for real nall=ncharged)
434
435
436 if(!fAODEvent && !fESDEvent)return;
437
438 //=================== AOD ===============
439
440 if(fAODEvent){//AOD loop
441
442 //reset global values
443 fNMcPrimAccept=0;// number of accepted primaries
444 fNRecAccept=0; // number of accepted tracks
445
446 if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec
447
448 if(!fMcOnly){
449 //step 5 = TrigVtxRecNrec
450 ntracks = LoopAOD(pt, eta, phi, charge, nTracksTracklets, 5);//read tracks
451 if(pt.size()){ //(internally ntracks=fNRecAccept)
452 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//analyse
453 }
454
455 if(fUseMC){
456 // step 4 = TrigVtxRecMcPropNrec
457 ntracks = LoopAODRecMcProp(pt, eta, phi, charge, nTracksTracklets, 4);//read tracks
458 if(pt.size()){//(internally ntracks=fNRecAccept)
459 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4);//analyse
460 }
461 }
9ed461df 462 }
9ed461df 463
7099c89a 464 if(fUseMC){
465 // step 3 = TrigVtxMcNrec
466 ntracks = LoopAODMC(pt, eta, phi, charge, nTracksTracklets, 3);//read tracks
467 if(pt.size()){//(internally ntracks=fNRecAccept)
468 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse
b9ad6f04 469 }
117f99e3 470 }
7099c89a 471
472
473 }//check event (true)
474
475 //reset values
476 fNMcPrimAccept=0;// number of accepted primaries
477 fNRecAccept=0; // number of accepted tracks
478
479 if(CheckEvent(false)){// all events, with and without reconstucted vertex
480 ntracks = LoopAOD (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events
481 ntracks = LoopAODMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks
482 if(pt.size()){
483 Analyse(pt, eta, phi, charge, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 2); // step 2 = TrigAllMcNrec
484
485 Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1); // step 1 = TrigAllMcNmc
117f99e3 486 }
7099c89a 487
488 //step 0 (include not triggered events) is not possible with AODs generated before October 2011
489 //step 0 can be implemented for new AODs
490
117f99e3 491 }
7099c89a 492 }//AOD loop
493
494
495 //=================== ESD ===============
496
497
498 if(fESDEvent){//ESD loop
499
500 //reset values
501 fNMcPrimAccept=0;// number of accepted primaries
502 fNRecAccept=0; // number of accepted tracks
503
504 // instead of task->SelectCollisionCandidate(mask) in AddTask macro
505 Bool_t isSelected = (((AliInputEventHandler*)
506 (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
507 ->IsEventSelected() & AliVEvent::kMB);
508
9ed461df 509
7099c89a 510 if(fDebug){
511 Printf("IsSelected = %d", isSelected);
512 Printf("CheckEvent(true)= %d", CheckEvent(true));
513 Printf("CheckEvent(false)= %d", CheckEvent(false));
514 }
515
516 //check trigger
517 if(isSelected){ // has offline trigger
518
519 if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec
520
521 if(!fMcOnly){
522 //step 5 = TrigVtxRecNrec
523 ntracks = LoopESD(pt, eta, phi, charge,nTracksTracklets, 5);//read tracks
524 if(pt.size()){ //(internally ntracks=fNRecAccept)
525 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//analyse
526 }
527
528 if(fUseMC){
529 // step 4 = TrigVtxRecMcPropNrec
530 ntracks = LoopESDRecMcProp(pt, eta, phi, charge, nTracksTracklets, 4);//read tracks
531 if(pt.size()){//(internally ntracks=fNRecAccept)
532 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4);//analyse
533 }
534 }
b9ad6f04 535 }
7099c89a 536
537 if(fUseMC){
538 // step 3 = TrigVtxMcNrec
539 ntracks = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 3);//read tracks
540 if(pt.size()){//(internally ntracks=fNRecAccept)
541 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse
542 }
543 }
544
545
546 }//check event (true)
547
548 //reset values
549 fNMcPrimAccept=0;// number of accepted primaries
550 fNRecAccept=0; // number of accepted tracks
551
552 if(CheckEvent(false)){// all events, with and without reconstucted vertex
553 ntracks = LoopESD (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events
554 ntracks = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks
555 if(pt.size()){
556 Analyse(pt, eta, phi, charge, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 2); // step 2 = TrigAllMcNrec
557
558 Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1); // step 1 = TrigAllMcNmc
559
560 Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0); //first part of step 0 // step 0 = AllAllMcNmc
561 }
562
563
564 }
565 }// triggered event
566
567 else { // not selected by physics selection task = not triggered
568 if(CheckEvent(false)){
569 ntracks = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 0);//read tracks
570 if(pt.size())Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0); //second part of step 0 // step 0 = AllAllMcNmc
117f99e3 571 }
117f99e3 572 }
7099c89a 573
574
575 }//ESD loop
576
577
578
579
117f99e3 580
581}
582
583
584//________________________________________________________________________
7099c89a 585Int_t AliAnalysisTaskMinijet::LoopESD( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
586 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
587 vector<Int_t> &nTracksTracklets, const Int_t step)
117f99e3 588{
589 // gives back the number of esd tracks and pointer to arrays with track
590 // properties (pt, eta, phi)
9ed461df 591 // Uses TPC tracks with SPD vertex from now on
117f99e3 592
7099c89a 593 ptArray.clear();
594 etaArray.clear();
595 phiArray.clear();
596 chargeArray.clear();
597 nTracksTracklets.clear();
598
599 const AliESDVertex* vtxSPD = fESDEvent->GetPrimaryVertexSPD(); // uses track or SPD vertexer
600 fVertexZ[step]->Fill(vtxSPD->GetZ());
601
117f99e3 602 // Retreive the number of all tracks for this event.
603 Int_t ntracks = fESDEvent->GetNumberOfTracks();
7099c89a 604 if(fDebug>1) Printf("all ESD tracks: %d", ntracks);
117f99e3 605
606 //first loop to check how many tracks are accepted
607 //------------------
608 Int_t nAcceptedTracks=0;
609 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
9ed461df 610
7099c89a 611 AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
9ed461df 612 if (!esdTrack) {
7099c89a 613 Error("LoopESD", "Could not receive track %d", iTracks);
117f99e3 614 continue;
615 }
7099c89a 616
9ed461df 617 // use TPC only tracks with non default SPD vertex
618 AliESDtrack *track = AliESDtrackCuts::
619 GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
620 if(!track) continue;
7099c89a 621 if(!fCuts->AcceptTrack(track)) {
622 delete track;
623 continue;// apply TPC track cuts before constrain to SPD vertex
624 }
9ed461df 625 if(track->Pt()>0.){
626 // only constrain tracks above threshold
627 AliExternalTrackParam exParam;
628 // take the B-field from the ESD, no 3D fieldMap available at this point
629 Bool_t relate = false;
630 relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
631 kVeryBig,&exParam);
632 if(!relate){
633 delete track;
634 continue;
635 }
636 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
637 exParam.GetCovariance());
638 }
7099c89a 639 else{
640 delete track;
641 continue;// only if tracks have pt<=0
642 }
9ed461df 643
7099c89a 644 if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax) {
645 ptArray.push_back(track->Pt());
646 etaArray.push_back(track->Eta());
647 phiArray.push_back(track->Phi());
648 chargeArray.push_back(track->Charge());
9ed461df 649 ++nAcceptedTracks;
7099c89a 650 fHistPt->Fill(track->Pt());
651 }
9ed461df 652
653 // TPC only track memory needs to be cleaned up
654 if(track)
655 delete track;
656
117f99e3 657 }
7099c89a 658
659 //need to be checked
660 if(nAcceptedTracks==0) return 0;
661
662 //tracklet loop
663 Int_t ntrackletsAccept=0;
664 AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
665 Int_t ntracklets = mult->GetNumberOfTracklets();
666 for(Int_t i=0;i< ntracklets;i++){
667 if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
668 ntrackletsAccept++;
669 }
670 }
671 nTracksTracklets.push_back(nAcceptedTracks);
672 nTracksTracklets.push_back(ntrackletsAccept);
673 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
674 //where here also neutral particles are counted.
9ed461df 675
676
7099c89a 677 fVzEvent=vtxSPD->GetZ(); // needed for correction map
678 if(step==5 || step ==2) fNRecAccept=nAcceptedTracks;
679 return fNRecAccept;
117f99e3 680
117f99e3 681
7099c89a 682}
9ed461df 683
7099c89a 684//________________________________________________________________________
685Int_t AliAnalysisTaskMinijet::LoopESDRecMcProp( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
686 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
687 vector<Int_t> &nTracksTracklets, const Int_t step)
688{
689 // gives back the number of esd tracks and pointer to arrays with track
690 // properties (pt, eta, phi) of mc particles if available
691 // Uses TPC tracks with SPD vertex from now on
9ed461df 692
7099c89a 693 ptArray.clear();
694 etaArray.clear();
695 phiArray.clear();
696 chargeArray.clear();
697 nTracksTracklets.clear();
9ed461df 698
7099c89a 699
700 AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
701 if (!mcEvent) {
702 Error("LoopESDRecMcProp", "Could not retrieve MC event");
703 return 0;
704 }
705 AliStack* stack = MCEvent()->Stack();
706 if(!stack) return 0;
707
708
709 // Retreive the number of all tracks for this event.
710 Int_t ntracks = fESDEvent->GetNumberOfTracks();
711 if(fDebug>1)Printf("all ESD tracks: %d", ntracks);
712
713 const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
714 fVertexZ[step]->Fill(vtxSPD->GetZ());
715
716 //track loop
717 Int_t nAcceptedTracks=0;
117f99e3 718 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
7099c89a 719
720 AliVParticle *vtrack = fESDEvent->GetTrack(iTracks);
9ed461df 721 AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
722 if (!esdTrack) {
7099c89a 723 Error("LoopESDRecMcProp", "Could not receive track %d", iTracks);
117f99e3 724 continue;
725 }
7099c89a 726
727 // use TPC only tracks with non default SPD vertex
9ed461df 728 AliESDtrack *track = AliESDtrackCuts::
729 GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
730 if(!track) continue;
7099c89a 731 if(!fCuts->AcceptTrack(track)) {
732 delete track;
733 continue;// apply TPC track cuts before constrain to SPD vertex
734 }
9ed461df 735 if(track->Pt()>0.){
736 // only constrain tracks above threshold
737 AliExternalTrackParam exParam;
738 // take the B-field from the ESD, no 3D fieldMap available at this point
739 Bool_t relate = false;
740 relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
741 kVeryBig,&exParam);
742 if(!relate){
743 delete track;
744 continue;
745 }
746 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
747 exParam.GetCovariance());
748 }
7099c89a 749 else{
750 delete track;
751 continue;// only if tracks have pt<=0
752 }
9ed461df 753
7099c89a 754 //count tracks, if available, use mc particle properties
755 if(vtrack->GetLabel()<0){
756 if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax){
757 ptArray.push_back(track->Pt());
758 etaArray.push_back(track->Eta());
759 phiArray.push_back(track->Phi());
760 chargeArray.push_back(track->Charge());
761 ++nAcceptedTracks;
762 }
117f99e3 763 }
7099c89a 764 else{
765 TParticle *partOfTrack = stack->Particle(vtrack->GetLabel());
766 if (TMath::Abs(partOfTrack->Eta())<fEtaCut && partOfTrack->Pt()>fPtMin && partOfTrack->Pt()<fPtMax) {
767 ptArray.push_back(partOfTrack->Pt());
768 etaArray.push_back(partOfTrack->Eta());
769 phiArray.push_back(partOfTrack->Phi());
770 chargeArray.push_back(vtrack->Charge());
771 ++nAcceptedTracks;
772 }
773 }
774
9ed461df 775 // TPC only track memory needs to be cleaned up
776 if(track)
777 delete track;
778
117f99e3 779 }
780
7099c89a 781 if(nAcceptedTracks==0) return 0;
117f99e3 782
783 //tracklet loop
784 Int_t ntrackletsAccept=0;
785 AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
786 Int_t ntracklets = mult->GetNumberOfTracklets();
787 for(Int_t i=0;i< ntracklets;i++){
7099c89a 788 if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
117f99e3 789 ntrackletsAccept++;
790 }
791 }
792
7099c89a 793 nTracksTracklets.push_back(nAcceptedTracks);
794 nTracksTracklets.push_back(ntrackletsAccept);
795 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
796 //where here also neutral particles are counted.
797
798
799 //get mc vertex for correction maps
800 AliGenEventHeader* header = MCEvent()->GenEventHeader();
801 TArrayF mcV;
802 header->PrimaryVertex(mcV);
803 fVzEvent= mcV[2];
804
805 return fNRecAccept; // give back reconstructed value
117f99e3 806
117f99e3 807
808}
809
7099c89a 810
811
812
117f99e3 813//________________________________________________________________________
7099c89a 814Int_t AliAnalysisTaskMinijet::LoopESDMC(vector<Float_t> &ptArray, vector<Float_t> &etaArray,
815 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
816 vector<Int_t> &nTracksTracklets, const Int_t step)
117f99e3 817{
818 // gives back the number of charged prim MC particle and pointer to arrays
819 // with particle properties (pt, eta, phi)
820
7099c89a 821 ptArray.clear();
822 etaArray.clear();
823 phiArray.clear();
824 chargeArray.clear();
825 nTracksTracklets.clear();
826
827 fNMcPrimAccept=0;
b9ad6f04 828
117f99e3 829 AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
830 if (!mcEvent) {
7099c89a 831 Error("LoopESDMC", "Could not retrieve MC event");
117f99e3 832 return 0;
833 }
834
9ed461df 835 AliStack* stack = MCEvent()->Stack();
836 if(!stack) return 0;
117f99e3 837
838 Int_t ntracks = mcEvent->GetNumberOfTracks();
7099c89a 839 if(fDebug>1)Printf("MC particles: %d", ntracks);
117f99e3 840
7099c89a 841 //vertex
842 AliGenEventHeader* header = MCEvent()->GenEventHeader();
843 TArrayF mcV;
844 header->PrimaryVertex(mcV);
845 Float_t vzMC = mcV[2];
846 fVertexZ[step]->Fill(vzMC);
117f99e3 847
848 //----------------------------------
849 //first track loop to check how many chared primary tracks are available
850 Int_t nChargedPrimaries=0;
b9ad6f04 851 Int_t nAllPrimaries=0;
852
853 Int_t nPseudoTracklets=0;
117f99e3 854 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
855 AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
856 if (!track) {
7099c89a 857 Error("LoopESDMC", "Could not receive track %d", iTracks);
117f99e3 858 continue;
859 }
b9ad6f04 860
861
862 if(//count also charged particles in case of fSelectParticles==2 (only neutral)
863 !SelectParticlePlusCharged(
7099c89a 864 track->Charge(),
865 track->PdgCode(),
866 stack->IsPhysicalPrimary(track->Label())
867 )
b9ad6f04 868 )
869 continue;
870
b9ad6f04 871 //count number of pseudo tracklets
7099c89a 872 if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.0) ++nPseudoTracklets; //0.035
117f99e3 873 //same cuts as on ESDtracks
7099c89a 874 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
b9ad6f04 875
876 //count all primaries
877 ++nAllPrimaries;
878 //count charged primaries
879 if (track->Charge() != 0) ++nChargedPrimaries;
880
7099c89a 881 if(fDebug>2) Printf("PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
b9ad6f04 882
883
117f99e3 884 }
885 //----------------------------------
886
7099c89a 887 if(fDebug>2){
888 Printf("All in acceptance=%d",nAllPrimaries);
889 Printf("Charged in acceptance =%d",nChargedPrimaries);
890 }
b9ad6f04 891
892 fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
117f99e3 893
b9ad6f04 894 if(nAllPrimaries==0) return 0;
117f99e3 895 if(nChargedPrimaries==0) return 0;
896
117f99e3 897 //track loop
7099c89a 898 //Int_t iChargedPiK=0;
117f99e3 899 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
900 AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
901 if (!track) {
7099c89a 902 Error("LoopESDMC", "Could not receive track %d", iTracks);
117f99e3 903 continue;
904 }
b9ad6f04 905
906 if(!SelectParticle(
907 track->Charge(),
908 track->PdgCode(),
909 stack->IsPhysicalPrimary(track->Label())
910 )
911 )
912 continue;
913
914
117f99e3 915 //same cuts as on ESDtracks
7099c89a 916 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
117f99e3 917
7099c89a 918 if(fDebug>2) Printf("After: PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
b9ad6f04 919
117f99e3 920
921 fHistPtMC->Fill(track->Pt());
922 //fills arrays with track properties
7099c89a 923 ptArray.push_back(track->Pt());
924 etaArray.push_back(track->Eta());
925 phiArray.push_back(track->Phi());
926 chargeArray.push_back(track->Charge());
927
928
117f99e3 929 } //track loop
930
7099c89a 931 nTracksTracklets.push_back(nChargedPrimaries);
932 nTracksTracklets.push_back(nPseudoTracklets);
b9ad6f04 933 if(fSelectParticles!=2){
7099c89a 934 nTracksTracklets.push_back(nAllPrimaries);
b9ad6f04 935 }
936 else{
7099c89a 937 nTracksTracklets.push_back(nAllPrimaries-nChargedPrimaries); // only neutral
b9ad6f04 938 }
9ed461df 939
7099c89a 940
9ed461df 941 fNMcPrimAccept=nChargedPrimaries;
7099c89a 942
943 if(step==1){
944 fNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
945 fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
946 }
947
948 fVzEvent= mcV[2];
949 return fNRecAccept;
117f99e3 950
951}
952
953//________________________________________________________________________
7099c89a 954Int_t AliAnalysisTaskMinijet::LoopAOD( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
955 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
956 vector<Int_t> &nTracksTracklets, const Int_t step)
117f99e3 957{
958 // gives back the number of AOD tracks and pointer to arrays with track
959 // properties (pt, eta, phi)
960
7099c89a 961 ptArray.clear();
962 etaArray.clear();
963 phiArray.clear();
964 chargeArray.clear();
965 nTracksTracklets.clear();
966
967 TClonesArray *mcArray=0x0;
968 if(fAnalysePrimOnly){
969 mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
970 }
971
972
973 AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
974 Double_t vzAOD=vertex->GetZ();
975 fVertexZ[step]->Fill(vzAOD);
117f99e3 976
977 // Retreive the number of tracks for this event.
978 Int_t ntracks = fAODEvent->GetNumberOfTracks();
7099c89a 979 if(fDebug>1) Printf("AOD tracks: %d", ntracks);
117f99e3 980
981
982 Int_t nAcceptedTracks=0;
983 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
984 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
985 if (!track) {
7099c89a 986 Error("LoopAOD", "Could not receive track %d", iTracks);
117f99e3 987 continue;
988 }
7099c89a 989
990 AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
991
992 //use only tracks from primaries
993 if(fAnalysePrimOnly){
994 if(vtrack->GetLabel()<0)continue;
995 if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
996 }
997
998 if(track->TestFilterBit(128) && TMath::Abs(track->Eta())<fEtaCut
999 && track->Pt()>fPtMin && track->Pt()<fPtMax){
1000
1001 nAcceptedTracks++;
1002
1003 // Printf("dca= %f", track->DCA());
1004 //save track properties in vector
1005 ptArray.push_back(track->Pt());
1006 etaArray.push_back(track->Eta());
1007 phiArray.push_back(track->Phi());
1008 chargeArray.push_back(track->Charge());
1009 fHistPt->Fill(track->Pt());
1010
1011 }
1012 }
1013 //need to check this option for MC
1014 if(nAcceptedTracks==0) return 0;
1015
1016
1017 //tracklet loop
1018 Int_t ntrackletsAccept=0;
1019 AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1020 for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1021 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut){
1022 ++ntrackletsAccept;
1023 }
117f99e3 1024 }
7099c89a 1025
1026
1027 nTracksTracklets.push_back(nAcceptedTracks);
1028 nTracksTracklets.push_back(ntrackletsAccept);
1029 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1030 //where here also neutral particles are counted.
1031
117f99e3 1032
7099c89a 1033 fVzEvent= vzAOD;
1034 if(step==5 || step==2)fNRecAccept = nAcceptedTracks; // needed for MC case //step5 = TrigVtxRecNrec
1035 return fNRecAccept; // at the moment, always return reconstructed multiplicity
117f99e3 1036
7099c89a 1037}
1038
1039//________________________________________________________________________
1040Int_t AliAnalysisTaskMinijet::LoopAODRecMcProp( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
1041 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1042 vector<Int_t> &nTracksTracklets, const Int_t step)
1043{
1044 // gives back the number of AOD tracks and pointer to arrays with track
1045 // properties (pt, eta, phi)
1046
1047
1048 ptArray.clear();
1049 etaArray.clear();
1050 phiArray.clear();
1051 chargeArray.clear();
1052 nTracksTracklets.clear();
b9ad6f04 1053
b9ad6f04 1054
7099c89a 1055 // Retreive the number of tracks for this event.
1056 Int_t ntracks = fAODEvent->GetNumberOfTracks();
1057 if(fDebug>1) Printf("AOD tracks: %d", ntracks);
b9ad6f04 1058
7099c89a 1059
1060 //get array of mc particles
1061 TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1062 FindListObject(AliAODMCParticle::StdBranchName());
1063 if(!mcArray){
1064 Printf("No MC particle branch found");
1065 return kFALSE;
1066 }
1067
1068 AliAODVertex* vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1069 Double_t vzAOD=vtx->GetZ();
1070 fVertexZ[step]->Fill(vzAOD);
117f99e3 1071
7099c89a 1072 Int_t nAcceptedTracks=0;
117f99e3 1073 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1074 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
7099c89a 1075
1076 AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1077
117f99e3 1078 if (!track) {
7099c89a 1079 Error("LoopAODRecMcProp", "Could not receive track %d", iTracks);
b9ad6f04 1080 continue;
117f99e3 1081 }
7099c89a 1082
1083 //use only tracks from primaries
1084 if(fAnalysePrimOnly){
1085 if(vtrack->GetLabel()<0)continue;
1086 if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1087 }
117f99e3 1088
7099c89a 1089 if(track->TestFilterBit(128) && TMath::Abs(track->Eta())<fEtaCut &&
1090 track->Pt()>fPtMin && track->Pt()<fPtMax){
1091
1092 nAcceptedTracks++;
1093
1094 //save track properties in vector
1095 if(vtrack->GetLabel()<0){ //fake tracks
1096 // Printf("Fake track");
1097 // continue;
1098 ptArray.push_back(track->Pt());
1099 etaArray.push_back(track->Eta());
1100 phiArray.push_back(track->Phi());
1101 chargeArray.push_back(track->Charge());
1102 }
1103 else{//mc properties
1104 AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1105
1106 ptArray.push_back(partOfTrack->Pt());
1107 etaArray.push_back(partOfTrack->Eta());
1108 phiArray.push_back(partOfTrack->Phi());
1109 chargeArray.push_back(vtrack->Charge());//partOfTrack?
1110 }
117f99e3 1111
7099c89a 1112 }
1113 }
1114 //need to check this option for MC
1115 if(nAcceptedTracks==0) return 0;
117f99e3 1116
1117 //tracklet loop
1118 Int_t ntrackletsAccept=0;
1119 AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1120 for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
7099c89a 1121 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut ){
117f99e3 1122 ++ntrackletsAccept;
1123 }
1124 }
1125
7099c89a 1126
1127 nTracksTracklets.push_back(nAcceptedTracks);
1128 nTracksTracklets.push_back(ntrackletsAccept);
1129 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1130 //where here also neutral particles are counted.
117f99e3 1131
117f99e3 1132
7099c89a 1133 //check vertex (mc)
1134 AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1135 FindListObject(AliAODMCHeader::StdBranchName());
1136 Float_t vzMC = aodMCheader->GetVtxZ();
1137
1138 fVzEvent= vzMC;
1139 return fNRecAccept;//this is the rec value from step 5
1140
1141}
1142
1143
117f99e3 1144
1145//________________________________________________________________________
7099c89a 1146Int_t AliAnalysisTaskMinijet::LoopAODMC( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
1147 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1148 vector<Int_t> &nTracksTracklets, const Int_t step)
117f99e3 1149{
b7c0438e 1150 // gives back the number of AOD MC particles and pointer to arrays with particle
1151 // properties (pt, eta, phi)
117f99e3 1152
7099c89a 1153 ptArray.clear();
1154 etaArray.clear();
1155 phiArray.clear();
1156 chargeArray.clear();
1157 nTracksTracklets.clear();
1158
1159 //check vertex
1160 AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1161 FindListObject(AliAODMCHeader::StdBranchName());
1162 Float_t vzMC = aodMCheader->GetVtxZ();
1163 fVertexZ[step]->Fill(vzMC);
1164
1165
b7c0438e 1166 //retreive MC particles from event
1167 TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1168 FindListObject(AliAODMCParticle::StdBranchName());
1169 if(!mcArray){
7099c89a 1170 Printf("No MC particle branch found");
b7c0438e 1171 return kFALSE;
117f99e3 1172 }
1173
b7c0438e 1174 Int_t ntracks = mcArray->GetEntriesFast();
7099c89a 1175 if(fDebug>1)Printf("MC particles: %d", ntracks);
117f99e3 1176
b7c0438e 1177
1178 // Track loop: chek how many particles will be accepted
7099c89a 1179 //Float_t vzMC=0.;
b9ad6f04 1180 Int_t nChargedPrim=0;
1181 Int_t nAllPrim=0;
1182 Int_t nPseudoTracklets=0;
b7c0438e 1183 for (Int_t it = 0; it < ntracks; it++) {
1184 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
117f99e3 1185 if (!track) {
7099c89a 1186 Error("LoopAODMC", "Could not receive particle %d", it);
117f99e3 1187 continue;
1188 }
b9ad6f04 1189
b18f7dfb 1190 if(!SelectParticlePlusCharged(
7099c89a 1191 track->Charge(),
1192 track->PdgCode(),
1193 track->IsPhysicalPrimary()
1194 )
b9ad6f04 1195 )
1196 continue;
b9ad6f04 1197
7099c89a 1198 if(TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.0)++nPseudoTracklets; //0.035
1199 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1200
b9ad6f04 1201 nAllPrim++;
1202 if(track->Charge()!=0) nChargedPrim++;
b18f7dfb 1203
b9ad6f04 1204 }
1205
b7c0438e 1206
b9ad6f04 1207 if(nAllPrim==0) return 0;
1208 if(nChargedPrim==0) return 0;
117f99e3 1209
7099c89a 1210 //generate array with size of number of accepted tracks
1211 fChargedPi0->Fill(nAllPrim,nChargedPrim);
1212
117f99e3 1213
b7c0438e 1214 // Track loop: fill arrays for accepted tracks
7099c89a 1215 // Int_t iChargedPrim=0;
b7c0438e 1216 for (Int_t it = 0; it < ntracks; it++) {
1217 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1218 if (!track) {
7099c89a 1219 Error("LoopAODMC", "Could not receive particle %d", it);
b7c0438e 1220 continue;
1221 }
b9ad6f04 1222
1223 if(!SelectParticle(
1224 track->Charge(),
1225 track->PdgCode(),
1226 track->IsPhysicalPrimary()
1227 )
1228 )
1229 continue;
1230
7099c89a 1231 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1232
b9ad6f04 1233 fHistPtMC->Fill(track->Pt());
7099c89a 1234 ptArray.push_back(track->Pt());
1235 etaArray.push_back(track->Eta());
1236 phiArray.push_back(track->Phi());
1237 chargeArray.push_back(track->Charge());
b7c0438e 1238 }
b9ad6f04 1239
7099c89a 1240 nTracksTracklets.push_back(nChargedPrim);
1241 nTracksTracklets.push_back(nPseudoTracklets);
b18f7dfb 1242 if(fSelectParticles!=2){
7099c89a 1243 nTracksTracklets.push_back(nAllPrim);
b18f7dfb 1244 }
1245 else{
7099c89a 1246 nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral
b18f7dfb 1247 }
7099c89a 1248
1249
b18f7dfb 1250
7099c89a 1251 fVzEvent= vzMC;
1252 fNMcPrimAccept = nChargedPrim;
1253 if(step==1){ // step 1 = Trig All Mc Nmc
1254 fNmcNch->Fill( fNMcPrimAccept,fNRecAccept);
1255 fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1256 }
1257 return fNRecAccept; // rec value from step 5 or step 2
1258
117f99e3 1259
1260}
1261
1262//________________________________________________________________________
7099c89a 1263void AliAnalysisTaskMinijet::Analyse(const vector<Float_t> &pt,
1264 const vector<Float_t> &eta,
1265 const vector<Float_t> &phi,
1266 const vector<Short_t> &charge,
1267 const Int_t ntracksCharged,
1268 const Int_t ntracklets,
1269 const Int_t nAll,
1270 const Int_t step)
117f99e3 1271{
1272
1273 // analyse track properties (comming from either ESDs or AODs) in order to compute
1274 // mini jet activity (chared tracks) as function of charged multiplicity
1275
7099c89a 1276 fStep->Fill(step);
117f99e3 1277
b18f7dfb 1278 if(fDebug){
7099c89a 1279 Printf("Analysis Step=%d", step);
1280 if(fDebug>2){
1281 Printf("nAll=%d",nAll);
1282 Printf("nCharged=%d",ntracksCharged);
1283 for (Int_t i = 0; i < nAll; i++) {
1284 Printf("pt[%d]=%f",i,pt[i]);
1285 }
1286 }
1287 }
1288
1289 //calculation of mean pt for all tracks in case of step==0
1290 if(step==5 || step ==2){ // rec step
1291 Double_t meanPt=0.;
1292 Double_t leadingPt=0.;
1293 for (Int_t i = 0; i < nAll; i++) {
1294 if(pt[i]<0.01)continue;
1295 meanPt+=pt[i];
1296 if(leadingPt<pt[i])leadingPt=pt[i];
1297 }
1298 meanPt=meanPt/nAll;
1299 fMeanPtRec=meanPt;
1300 fLeadingPtRec=leadingPt;
b18f7dfb 1301 }
7099c89a 1302
1303 Double_t propEvent[4] = {ntracksCharged,fVzEvent,fMeanPtRec, fLeadingPtRec}; //vz: {rec, mc, mc}, meanPt and Nrec is always rec value
1304 fMapEvent[step]->Fill(propEvent);
1305
1306 fNcharge[step]->Fill(ntracksCharged);
117f99e3 1307
7099c89a 1308 Float_t ptEventAxis=0; // pt event axis
1309 Float_t etaEventAxis=0; // eta event axis
1310 Float_t phiEventAxis=0; // phi event axis
1311 Short_t chargeEventAxis=0; // charge event axis
117f99e3 1312
1313 Float_t ptOthers = 0; // pt others // for all other tracks around event axis -> see loop
1314 Float_t etaOthers = 0; // eta others
1315 Float_t phiOthers = 0; // phi others
7099c89a 1316 Short_t chargeOthers = 0; // charge others
b18f7dfb 1317
a640cfb3 1318 Int_t *pindexInnerEta = new Int_t [nAll+1];
1319 Float_t *ptInnerEta = new Float_t[nAll+1];
b18f7dfb 1320
7099c89a 1321
1322
b9ad6f04 1323 for (Int_t i = 0; i < nAll; i++) {
7099c89a 1324
1325 if(pt[i]<0.01)continue;
1326
1327 //fill single particle correction for first step of pair correction
1328 Double_t propAll[3] = {eta[i],pt[i],ntracksCharged };
1329 fMapAll[step]->Fill(propAll);
1330
1331
117f99e3 1332 //filling of simple check plots
7099c89a 1333 if(pt[i]<0.7) continue;
1334 fPt[step] -> Fill( pt[i]);
1335 fEta[step] -> Fill(eta[i]);
1336 fPhi[step] -> Fill(phi[i]);
1337 fPhiEta[step]-> Fill(phi[i], eta[i]);
9ed461df 1338
b9ad6f04 1339 pindexInnerEta[i]=0; //set all values to zero
1340 //fill new array for eta check
1341 ptInnerEta[i]= pt[i];
1342
117f99e3 1343 }
7099c89a 1344
1345
b18f7dfb 1346
117f99e3 1347 // define event axis: leading or random track (with pt>fTriggerPtCut)
1348 // ---------------------------------------
1349 Int_t highPtTracks=0;
b9ad6f04 1350 Int_t highPtTracksInnerEta=0;
7099c89a 1351 // Double_t highPtTracksInnerEtaCorr=0;
9ed461df 1352 Int_t mult09=0;
b9ad6f04 1353
1354 //count high pt tracks and high pt tracks in acceptance for seeds
1355 for(Int_t i=0;i<nAll;i++){
1356
7099c89a 1357 if(pt[i]<0.01)continue;
1358
9ed461df 1359 if(TMath::Abs(eta[i])<0.9){
1360 mult09++;
1361 }
1362
b9ad6f04 1363 if(pt[i]>fTriggerPtCut) {
1364 highPtTracks++;
1365 }
1366
1367 // seed should be place in middle of acceptance, that complete cone is in acceptance
9ed461df 1368 if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
b9ad6f04 1369
1370 // Printf("eta=%f", eta[i]);
1371 highPtTracksInnerEta++;
7099c89a 1372
b9ad6f04 1373 }
1374 else{
1375 ptInnerEta[i]=0;
1376 }
117f99e3 1377 }
1378
b9ad6f04 1379
1380 //sort array in order to get highest pt tracks first
1381 //index can be computed with pindexInnerEta[number]
1382 if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1383
7099c89a 1384 // plot of multiplicity distributions
1385 fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1386 fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1387
117f99e3 1388 if(ntracklets){
7099c89a 1389 fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1390 fNchTracklet[step]->Fill(ntracklets, ntracksCharged);
1391 fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
117f99e3 1392 }
1393
1394 //analysis can only be performed with event axis, defined by high pt track
7099c89a 1395
b9ad6f04 1396
1397 if(highPtTracks>0 && highPtTracksInnerEta>0){
1398
7099c89a 1399 // build pairs in two track loops
9ed461df 1400 // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
7099c89a 1401 for(Int_t axis=0;(axis<nAll) && (pt[pindexInnerEta[axis]]>=fTriggerPtCut); axis++){
1402
117f99e3 1403 //EventAxisRandom track properties
b9ad6f04 1404 ptEventAxis = pt [pindexInnerEta[axis]];
1405 etaEventAxis = eta[pindexInnerEta[axis]];
1406 phiEventAxis = phi[pindexInnerEta[axis]];
7099c89a 1407 chargeEventAxis = charge[pindexInnerEta[axis]];
1408 fPtSeed[step] -> Fill( ptEventAxis);
1409 fEtaSeed[step] -> Fill(etaEventAxis);
1410 fPhiSeed[step] -> Fill(phiEventAxis);
117f99e3 1411
117f99e3 1412
7099c89a 1413 Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged };
1414 fMapSingleTrig[step]->Fill(prop);
9ed461df 1415
7099c89a 1416 //associated tracks
1417 for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) {
117f99e3 1418
7099c89a 1419 if(pt[pindexInnerEta[iTrack]]<fAssociatePtCut) continue;
9ed461df 1420
7099c89a 1421 if(fSelectParticlesAssoc==1){
1422 if(charge[pindexInnerEta[iTrack]]==0)continue;
1423 }
1424 if(fSelectParticlesAssoc==2){
1425 if(charge[pindexInnerEta[iTrack]]!=0)continue;
1426 }
9ed461df 1427
9ed461df 1428
7099c89a 1429 ptOthers = pt [pindexInnerEta[iTrack]];
1430 etaOthers = eta[pindexInnerEta[iTrack]];
1431 phiOthers = phi[pindexInnerEta[iTrack]];
1432 chargeOthers = charge[pindexInnerEta[iTrack]];
1433
1434
1435 //plot only properties of tracks with pt>ptassoc
1436 fPtOthers[step] -> Fill( ptOthers);
1437 fEtaOthers[step] -> Fill(etaOthers);
1438 fPhiOthers[step] -> Fill(phiOthers);
1439 fPtEtaOthers[step] -> Fill(ptOthers, etaOthers);
117f99e3 1440
7099c89a 1441 // if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]);
9ed461df 1442
7099c89a 1443 Float_t dPhi = phiOthers-phiEventAxis;
1444 if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1445 else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1446 Float_t dEta=etaOthers-etaEventAxis;
9ed461df 1447
1448
7099c89a 1449 fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta);
1450 fDPhiEventAxis[step]->Fill(dPhi);
1451
1452 //check outliers
1453 if(ptEventAxis< 0.4 || ptEventAxis > 100) Printf("particles out of range pt");
1454 if(ntracksCharged<0 || ntracksCharged>150) Printf("particles out of range ncharge");
1455 if(TMath::Abs(dEta)>2*fEtaCut) {
1456 Printf("particles out of range dEta");
1457 Printf("eta1=%f, eta2=%f", etaOthers, etaEventAxis);
1458 Printf("step=%d",step);
1459 }
1460 if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){
1461 Printf("particles out of range dPhi");
1462 Printf("phi1=%f, phi2=%f", phiOthers, phiEventAxis);
1463 }
1464
1465 Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers);
1466
1467 Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, isLikeSign };
1468 fMapPair[step]->Fill(prop6);
1469
1470 }// end of inner track loop
1471
1472 } //end of outer track loop
1473
1474 fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1475 fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1476 fTriggerTracklet[step]->Fill(ntracklets);
117f99e3 1477
9ed461df 1478
117f99e3 1479 }//if there is at least one high pt track
1480
b9ad6f04 1481
1482 if(pindexInnerEta){// clean up array memory used for TMath::Sort
1483 delete[] pindexInnerEta;
1484 pindexInnerEta=0;
117f99e3 1485 }
1486
b9ad6f04 1487 if(ptInnerEta){// clean up array memory used for TMath::Sort
1488 delete[] ptInnerEta;
1489 ptInnerEta=0;
1490 }
117f99e3 1491
1492}
1493
7099c89a 1494
1495
117f99e3 1496//________________________________________________________________________
1497void AliAnalysisTaskMinijet::Terminate(Option_t*)
1498{
1499 //terminate function is called at the end
1500 //can be used to draw histograms etc.
1501
7099c89a 1502
117f99e3 1503}
1504
1505//________________________________________________________________________
a640cfb3 1506Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, Bool_t prim)
b9ad6f04 1507{
1508 //selection of mc particle
1509 //fSelectParticles=0: use charged primaries and pi0 and k0
1510 //fSelectParticles=1: use only charged primaries
1511 //fSelectParticles=2: use only pi0 and k0
1512
1513 if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1514 if(charge==0){
1515 if(!(pdg==111||pdg==130||pdg==310))
1516 return false;
1517 }
1518 else{// charge !=0
1519 if(!prim)
1520 return false;
1521 }
1522 }
1523
1524 else if(fSelectParticles==1){
1525 if (charge==0 || !prim){
1526 return false;
1527 }
1528 }
1529
1530 else{
1531 Printf("Error: wrong selection of charged/pi0/k0");
1532 return 0;
1533 }
1534
1535 return true;
1536
1537}
1538
1539//________________________________________________________________________
a640cfb3 1540Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
b9ad6f04 1541{
1542 //selection of mc particle
1543 //fSelectParticles=0: use charged primaries and pi0 and k0
1544 //fSelectParticles=1: use only charged primaries
1545 //fSelectParticles=2: use only pi0 and k0
1546
1547 if(fSelectParticles==0){
1548 if(charge==0){
1549 if(!(pdg==111||pdg==130||pdg==310))
1550 return false;
1551 }
1552 else{// charge !=0
1553 if(!prim)
1554 return false;
1555 }
7099c89a 1556 }
b9ad6f04 1557 else if (fSelectParticles==1){
1558 if (charge==0 || !prim){
1559 return false;
1560 }
1561 }
1562 else if(fSelectParticles==2){
1563 if(!(pdg==111||pdg==130||pdg==310))
1564 return false;
1565 }
1566
1567 return true;
1568
1569}
1570
1571//________________________________________________________________________
a640cfb3 1572Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex)
117f99e3 1573{
7099c89a 1574 // This function tests the quality of an event (ESD/AOD) (rec and/or mc part)
1575 // recVertex=false: check if Mc events and stack is there, Nmc>0
1576 // recVertex=false: " + check if there is a good, reconstructed SPD vertex
1577 // defined by |z|<fVertexCut(10cm), Contributer>0, no PileUpFromSPD(3,0,8)
1578
1579
1580 if(fMode==0){//esd
1581
1582 //mc
1583 if(fUseMC){
1584
1585 //mc event
1586 AliMCEvent *mcEvente = (AliMCEvent*) MCEvent();
1587 if (!mcEvente) {
1588 Error("CheckEvent", "Could not retrieve MC event");
1589 return false;
1590 }
1591
1592 //stack
1593 AliStack* stackg = MCEvent()->Stack();
1594 if(!stackg) return false;
1595 Int_t ntracksg = mcEvente->GetNumberOfTracks();
1596 if(ntracksg<0) return false;
1597
1598 //vertex
1599 //AliGenEventHeader* headerg = MCEvent()->GenEventHeader();
1600 //TArrayF mcVg;
1601 //headerg->PrimaryVertex(mcVg);
1602 // if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 &&
1603 // TMath::Abs(mcVg[0])<1e-8) return false;
1604 // Float_t vzMCg = mcVg[2];
1605 // if(TMath::Abs(vzMCg)>fVertexZCut) return false;
1606 }
1607
1608 //rec
1609 if(recVertex==true){
1610 if(fESDEvent->IsPileupFromSPD(3,0,8))return false;
1611
1612 //rec vertex
1613 // const AliESDVertex* vertexESDg = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
1614 // if(!vertexESDg) return false;
1615 // if(vertexESDg->GetNContributors()<=0)return false;
1616 // Float_t fVzg= vertexESDg->GetZ();
1617 // if(TMath::Abs(fVzg)>fVertexZCut) return false;
1618
1619 //rec spd vertex
1620 const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
1621 if(!vtxSPD) return false;
1622 if(vtxSPD->GetNContributors()<=0)return false;
1623 Float_t fVzSPD= vtxSPD->GetZ();
1624 if(TMath::Abs(fVzSPD)>fVertexZCut) return false;
1625 }
1626 return true;
117f99e3 1627
117f99e3 1628 }
7099c89a 1629
1630
1631 else if(fMode==1){ //aod
1632
1633 if(fUseMC){
1634 //mc
1635 // AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1636 // FindListObject(AliAODMCHeader::StdBranchName());
1637 // Float_t vzMC = aodMCheader->GetVtxZ();
1638 // if(TMath::Abs(vzMC)>fVertexZCut) return false;
1639
1640 //retreive MC particles from event
1641 TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1642 FindListObject(AliAODMCParticle::StdBranchName());
1643 if(!mcArray){
1644 Printf("No MC particle branch found");
1645 return false;
1646 }
1647 }
1648
1649 //rec
1650 if(recVertex==true){
1651 if(fAODEvent->IsPileupFromSPD(3,0.8))return false;
1652
1653 // AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
1654 // if(!vertex) return false;
1655 // if(vertex->GetNContributors()<=0) return false;
1656 // Double_t vzAOD=vertex->GetZ();
1657 // if(TMath::Abs(vzAOD)<1e-9) return false;
1658 // if(TMath::Abs(vzAOD)>fVertexZCut) return false;
1659
1660 AliAODVertex* vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();
1661 if(!vertexSPD) return false;
1662 if(vertexSPD->GetNContributors()<=0) return false;
1663 Double_t vzSPD=vertexSPD->GetZ();
1664 //if(TMath::Abs(vzSPD)<1e-9) return false;
1665 if(TMath::Abs(vzSPD)>fVertexZCut) return false;
1666 }
1667 return true;
1668
117f99e3 1669 }
7099c89a 1670
1671 else {
1672 Printf("Wrong mode!");
1673 return false;
9ed461df 1674 }
7099c89a 1675
1676}
1677
1678//_____________________________________________________________________________
1679const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins,
1680 const Double_t xmin,
1681 const Double_t xmax)
1682{
1683 // returns pointer to an array with can be used to build a logarithmic axis
1684 // it is user responsibility to delete the array
1685
1686 Double_t logxmin = TMath::Log10(xmin);
1687 Double_t logxmax = TMath::Log10(xmax);
1688 Double_t binwidth = (logxmax-logxmin)/nbins;
1689
1690 Double_t *xbins = new Double_t[nbins+1];
1691
1692 xbins[0] = xmin;
1693 for (Int_t i=1;i<=nbins;i++) {
1694 xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth);
117f99e3 1695 }
1696
7099c89a 1697 return xbins;
117f99e3 1698}
b9ad6f04 1699
7099c89a 1700//_____________________________________________________________________________
a640cfb3 1701Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis,
1702 const Short_t chargeOthers)
7099c89a 1703{
1704 // compute if charge of two particles/tracks has same sign or different sign
1705
1706 if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers);
1707
1708 if(chargeEventAxis<0){
1709 if(chargeOthers<0){
1710 return true;
1711 }
1712 else if(chargeOthers>0){
1713 return false;
1714 }
1715 }
1716
1717 else if(chargeEventAxis>0){
1718 if(chargeOthers>0){
1719 return true;
1720 }
1721 else if(chargeOthers<0){
1722 return false;
1723 }
1724 }
1725
1726 else{
1727 Printf("Error: Charge not lower nor higher as zero");
1728 return false;
1729 }
1730
1731 Printf("Error: Check values of Charge");
1732 return false;
1733}
1734