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