12 #include "AliVEvent.h"
13 #include "AliVParticle.h"
15 #include "AliESDEvent.h"
16 #include "AliESDtrack.h"
17 #include "AliMultiplicity.h"
18 #include "AliESDtrackCuts.h"
20 #include "AliAODEvent.h"
21 #include "AliAODTrack.h"
22 #include "AliAODMCParticle.h"
23 #include "AliAODMCHeader.h"
26 #include "AliMCEvent.h"
27 #include "AliMCParticle.h"
28 #include "AliGenEventHeader.h"
29 #include "AliCentrality.h"
31 #include "AliAnalysisManager.h"
32 #include "AliInputEventHandler.h"
39 #include "AliAnalysisTaskMinijet.h"
41 // Analysis task for two-particle correlations using all particles over pt threshold
42 // pt_trig threshold for trigger particle (event axis) and pt_assoc for possible associated particles.
43 // Extract mini-jet yield and fragmentation properties via Delta-Phi histograms of these correlations
44 // post processing of analysis output via macro plot3and2Gaus.C
45 // Can use ESD or AOD, reconstructed and Monte Carlo data as input
46 // Author: Eva Sicking, modifications by Emilia Leogrande
49 ClassImp(AliAnalysisTaskMinijet)
51 //________________________________________________________________________
52 AliAnalysisTaskMinijet::AliAnalysisTaskMinijet(const char *name)
53 : AliAnalysisTaskSE(name),
57 fAnalysePrimOnly(kFALSE),// not used
70 fSelectParticlesAssoc(1),
73 fCorrStrangeness(true),
74 fThreeParticleCorr(false),
81 fNRecAcceptStrangeCorr(0),
82 fNMcPrimAcceptTracklet(0),
83 fNRecAcceptTracklet(0),
95 fCorruptedChunksAfter(0),
99 fNmcNchVtxStrangeCorr(0),
103 fNmcNchVtxTracklet(0),
104 fPNmcNchVtxTracklet(0),
108 fCentralityMethod("")
113 for(Int_t i = 0;i< 8;i++){
114 fMapSingleTrig[i] = 0;
140 fDPhiDEtaEventAxis[i] = 0;
141 fDPhiDEtaEventAxisSeeds[i]= 0;
144 fTriggerNchSeeds[i] = 0;
145 fTriggerTracklet[i] = 0;
150 fNch07Tracklet[i] = 0;
152 fPNch07Tracklet[i] = 0;
153 fDPhiEventAxis[i] = 0;
156 DefineOutput(1, TList::Class());
159 //________________________________________________________________________
160 AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet()
164 if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists;
167 //________________________________________________________________________
168 void AliAnalysisTaskMinijet::UserCreateOutputObjects()
172 if(fDebug) Printf("In User Create Output Objects.");
174 Int_t nbinsCentr = 101;
175 Float_t minbinCentr=-0.5, maxbinCentr=100.5;
177 fStep = new TH1F("fStep", "fStep", 10, -0.5, 9.5);
178 fEventStat = new TH1F("fEventStat", "fEventStat", 10, -0.5, 9.5);
179 fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 150, 0.1, 3.1);
180 fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
181 fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
182 fHistPt->SetMarkerStyle(kFullCircle);
183 fNContrNtracklets = new TH2F ("fNContrNtracklets", ";N_{tracklets};N_{vtx contrib}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
184 fNContrNtracks = new TH2F ("fNContrNtracks", ";N_{tracks};N_{vtx contrib}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
185 fCorruptedChunks = new TH2F ("fCorruptedChunks",
186 ";N_{tracks,TPC};N_{tracks,ITS-TPC}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
187 fCorruptedChunksAfter = new TH2F ("fCorruptedChunksAfter",
188 ";N_{tracks,TPC};N_{tracks,ITS-TPC}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
194 fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 150, 0.1, 3.1);
195 fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)");
196 fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
197 fHistPtMC->SetMarkerStyle(kFullCircle);
199 fNmcNch = new TH2F("fNmcNch", "fNmcNch", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
200 fPNmcNch = new TProfile("pNmcNch", "pNmcNch", nbinsCentr,minbinCentr,maxbinCentr);
201 fNmcNchVtx = new TH2F("fNmcNchVtx", "fNmcNchVtx", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
202 fNmcNchVtxStrangeCorr = new TH2F("fNmcNchVtxStrangeCorr", "fNmcNchVtxStrangeCorr", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
203 fPNmcNchVtx = new TProfile("pNmcNchVtx", "pNmcNchVtx", nbinsCentr,minbinCentr,maxbinCentr);
205 fNmcNchTracklet = new TH2F("fNmcNchTracklet", "fNmcNchTracklet", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
206 fPNmcNchTracklet = new TProfile("pNmcNchTracklet", "pNmcNchTracklet", nbinsCentr,minbinCentr,maxbinCentr);
207 fNmcNchVtxTracklet = new TH2F("fNmcNchVtxTracklet", "fNmcNchVtxTracklet", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
208 fPNmcNchVtxTracklet = new TProfile("pNmcNchVtxTracklet", "pNmcNchVtxTracklet", nbinsCentr,minbinCentr,maxbinCentr);
214 fChargedPi0 = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5);
215 fVertexCheck = new TH1F("fVertexCheck", "fVertexCheck", 100, -0.5, 99.5);
216 fPropagateDca = new TH1F("fPropagateDca", "fPropagateDca", 2, -0.5, 1.5);
218 //----------------------
219 //bins for pt in THnSpare
220 Double_t ptMin = 0.0, ptMax = 100.;
222 Double_t binsPt[] = {0.0, 0.1, 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8,
223 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,
224 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};
226 //Int_t nPtBinsAll = 100;
227 //Double_t ptMinAll = 1.e-2, ptMaxAll = 100.;
228 //Double_t *binsPtAll = 0;
229 //binsPtAll = (Double_t*)CreateLogAxis(nPtBinsAll,ptMinAll,ptMaxAll);
232 // Double_t ptMin2 = 0.0, ptMax2 = 100.;
233 // Int_t nPtBins2 = 9;
234 // Double_t binsPt2[] = {0.1, 0.4, 0.7, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0};
236 // Int_t nPtBins2 = 10;
237 // Double_t ptMin2 = 0.4, ptMax2 = 100.;
238 // Double_t *binsPt2 = 0;
239 // binsPt2 = CreateLogAxis(nPtBins2,ptMin2,ptMax2);
242 Int_t binsEffHisto[3] = {Int_t(fEtaCut*20), nPtBins, nbinsCentr };
243 Double_t minEffHisto[3] = {-fEtaCut, ptMin, minbinCentr };
244 Double_t maxEffHisto[3] = {fEtaCut, ptMax, maxbinCentr};
247 Int_t binsEffHisto5[6] = { nPtBins, nPtBins, 36, 90, nbinsCentr , 2 };
248 Double_t minEffHisto5[6] = { ptMin, ptMin, -2*fEtaCut, -0.5*TMath::Pi(), minbinCentr , -0.5 };
249 Double_t maxEffHisto5[6] = { ptMax, ptMax, 2*fEtaCut, 1.5*TMath::Pi(), maxbinCentr , 1.5 };
253 Int_t binsEvent[4] = { nbinsCentr, 20, 50, nPtBins };
254 Double_t minEvent[4] = { minbinCentr, -10, 0, ptMin };
255 Double_t maxEvent[4] = { maxbinCentr, 10, 10, ptMax };
258 Int_t binsAll[3] = {Int_t(fEtaCut*20), nPtBins, nbinsCentr };
259 Double_t minAll[3] = {-fEtaCut, ptMin, minbinCentr };
260 Double_t maxAll[3] = {fEtaCut, ptMax, maxbinCentr };
263 Int_t binsThree[3] = { 90, 90, nbinsCentr};
264 Double_t minThree[3] = {-0.5*TMath::Pi(), -0.5*TMath::Pi(), minbinCentr};
265 Double_t maxThree[3] = { 1.5*TMath::Pi(), 1.5*TMath::Pi(), maxbinCentr};
271 //--------------------
272 TString dataType[2] ={"ESD", "AOD"};
273 TString labels[8]={Form("%sAllAllMcNmc",dataType[fMode].Data()),
274 Form("%sTrigAllMcNmc",dataType[fMode].Data()),
275 Form("%sTrigVtxMcNmc",dataType[fMode].Data()),
276 Form("%sTrigVtxMcNrec",dataType[fMode].Data()),
277 Form("%sTrigVtxRecMcPropNrec",dataType[fMode].Data()),
278 Form("%sTrigVtxRecNrec",dataType[fMode].Data()),
279 Form("%sTrigVtxRecMcPropNrecStrangeCorr",dataType[fMode].Data()),
280 Form("%sTrigVtxRecNrecStrangeCorr",dataType[fMode].Data())};
283 for(Int_t i=0;i<8;i++){
285 fMapSingleTrig[i] = new THnSparseD(Form("fMapSingleTrig%s", labels[i].Data()),"eta:pt:Nrec",
286 3,binsEffHisto,minEffHisto,maxEffHisto);
287 fMapSingleTrig[i]->SetBinEdges(1,binsPt);
288 fMapSingleTrig[i]->GetAxis(0)->SetTitle("#eta");
289 fMapSingleTrig[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
290 fMapSingleTrig[i]->GetAxis(2)->SetTitle("N_{rec}");
291 fMapSingleTrig[i]->Sumw2();
293 fMapPair[i] = new THnSparseD(Form("fMapPair%s", labels[i].Data()),"pt_trig:pt_assoc:DeltaEta:DeltaPhi:Nrec:likesign",
294 6,binsEffHisto5,minEffHisto5,maxEffHisto5);
295 fMapPair[i]->SetBinEdges(0,binsPt);
296 fMapPair[i]->SetBinEdges(1,binsPt);
297 fMapPair[i]->GetAxis(0)->SetTitle("p_{T, trig} (GeV/c)");
298 fMapPair[i]->GetAxis(1)->SetTitle("p_{T, assoc} (GeV/c)");
299 fMapPair[i]->GetAxis(2)->SetTitle("#Delta #eta");
300 fMapPair[i]->GetAxis(3)->SetTitle("#Delta #phi");
301 fMapPair[i]->GetAxis(4)->SetTitle("N_{rec}");
302 fMapPair[i]->GetAxis(5)->SetTitle("Like-sign or Unlike-sign");
303 fMapPair[i]->Sumw2();
306 fMapEvent[i] = new THnSparseD(Form("fMapEvent%s", labels[i].Data()),"Nrec:vertexZ:meanPt:leadingPt",
307 4,binsEvent,minEvent,maxEvent);
308 fMapEvent[i]->GetAxis(0)->SetTitle("N_{rec}");
309 fMapEvent[i]->GetAxis(1)->SetTitle("z_{vertex} (cm)");
310 fMapEvent[i]->GetAxis(2)->SetTitle("meanPt");
311 fMapEvent[i]->SetBinEdges(3,binsPt);
312 fMapEvent[i]->GetAxis(3)->SetTitle("leadingPt");
313 fMapEvent[i]->Sumw2();
315 fMapAll[i] = new THnSparseD(Form("fMapAll%s", labels[i].Data()),"eta:pt:Nrec",
316 3,binsAll,minAll,maxAll);
317 fMapAll[i]->SetBinEdges(1,binsPt);
318 fMapAll[i]->GetAxis(0)->SetTitle("#eta");
319 fMapAll[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
320 fMapAll[i]->GetAxis(2)->SetTitle("N_{rec}");
324 fMapThree[i] = new THnSparseD(Form("fMapThree%s", labels[i].Data()),"dphi1:dphi2:Nrec",
325 3,binsThree,minThree,maxThree);
326 fMapThree[i]->GetAxis(0)->SetTitle("#Delta#varphi_{1}");
327 fMapThree[i]->GetAxis(1)->SetTitle("#Delta#varphi_{2}");
328 fMapThree[i]->GetAxis(2)->SetTitle("N_{rec}");
329 fMapThree[i]->Sumw2();
332 fVertexZ[i] = new TH1F(Form("fVertexZ%s",labels[i].Data()),
333 Form("fVertexZ%s",labels[i].Data()) ,
335 fPt[i] = new TH1F(Form("fPt%s",labels[i].Data()),
336 Form("fPt%s",labels[i].Data()) ,
338 fNcharge[i] = new TH1F(Form("fNcharge%s",labels[i].Data()),
339 Form("fNcharge%s",labels[i].Data()) ,
341 fEta[i] = new TH1F (Form("fEta%s",labels[i].Data()),
342 Form("fEta%s",labels[i].Data()) ,
344 fPhi[i] = new TH1F(Form("fPhi%s",labels[i].Data()),
345 Form("fPhi%s",labels[i].Data()) ,
346 360, 0.,2*TMath::Pi());
347 fDcaXY[i] = new TH1F(Form("fDcaXY%s",labels[i].Data()),
348 Form("fDcaXY%s",labels[i].Data()) ,
350 fDcaZ[i] = new TH1F(Form("fDcaZ%s",labels[i].Data()),
351 Form("fDcaZ%s",labels[i].Data()) ,
353 fPtSeed[i] = new TH1F(Form("fPSeedt%s",labels[i].Data()),
354 Form("fPtSeed%s",labels[i].Data()) ,
356 fEtaSeed[i] = new TH1F (Form("fEtaSeed%s",labels[i].Data()),
357 Form("fEtaSeed%s",labels[i].Data()) ,
359 fPhiSeed[i] = new TH1F(Form("fPhiSeed%s",labels[i].Data()),
360 Form("fPhiSeed%s",labels[i].Data()) ,
361 360, 0.,2*TMath::Pi());
363 fPtOthers[i] = new TH1F(Form("fPOtherst%s",labels[i].Data()),
364 Form("fPtOthers%s",labels[i].Data()) ,
366 fEtaOthers[i] = new TH1F (Form("fEtaOthers%s",labels[i].Data()),
367 Form("fEtaOthers%s",labels[i].Data()) ,
369 fPhiOthers[i] = new TH1F(Form("fPhiOthers%s",labels[i].Data()),
370 Form("fPhiOthers%s",labels[i].Data()) ,
371 360, 0.,2*TMath::Pi());
372 fPtEtaOthers[i] = new TH2F(Form("fPtEtaOthers%s",labels[i].Data()),
373 Form("fPtEtaOthers%s",labels[i].Data()) ,
374 500, 0., 50, 100, -1., 1);
375 fPhiEta[i] = new TH2F(Form("fPhiEta%s",labels[i].Data()),
376 Form("fPhiEta%s",labels[i].Data()) ,
377 180, 0., 2*TMath::Pi(), 100, -1.,1.);
378 fDPhiDEtaEventAxis[i] = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()),
379 Form("fDPhiDEtaEventAxis%s",labels[i].Data()) ,
380 180, -0.5* TMath::Pi(), 1.5*TMath::Pi(), 9, -1.8,1.8);
381 fTriggerNch[i] = new TH1F(Form("fTriggerNch%s",labels[i].Data()),
382 Form("fTriggerNch%s",labels[i].Data()) ,
384 fTriggerNchSeeds[i] = new TH2F(Form("fTriggerNchSeeds%s",labels[i].Data()),
385 Form("fTriggerNchSeeds%s",labels[i].Data()) ,
386 250, -0.5, 249.5, 100, -0.5, 99.5);
387 fTriggerTracklet[i] = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()),
388 Form("fTriggerTracklet%s",labels[i].Data()) ,
390 fNch07Nch[i] = new TH2F(Form("fNch07Nch%s",labels[i].Data()),
391 Form("fNch07Nch%s",labels[i].Data()) ,
392 250, -2.5, 247.5,250, -2.5, 247.5);
393 fPNch07Nch[i] = new TProfile(Form("pNch07Nch%s",labels[i].Data()),
394 Form("pNch07Nch%s",labels[i].Data()) ,
396 fNch07Tracklet[i] = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()),
397 Form("fNch07Tracklet%s",labels[i].Data()) ,
398 250, -2.5, 247.5,250, -2.5, 247.5);
399 fNchTracklet[i] = new TH2F(Form("fNchTracklet%s",labels[i].Data()),
400 Form("fNchTracklet%s",labels[i].Data()) ,
401 250, -2.5, 247.5,250, -2.5, 247.5);
402 fPNch07Tracklet[i] = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()),
403 Form("pNch07Tracklet%s",labels[i].Data()) ,
405 fDPhiEventAxis[i] = new TH1F(Form("fDPhiEventAxis%s",
407 Form("fDPhiEventAxis%s",
409 180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
410 fDPhi1DPhi2[i] = new TH2F(Form("fDPhi1DPhi2%s",labels[i].Data()),
411 Form("fDPhi1DPhi2%s",labels[i].Data()),
412 180, -0.5* TMath::Pi(), 1.5*TMath::Pi(),
413 180, -0.5* TMath::Pi(), 1.5*TMath::Pi());
416 fHists = new TList();
420 fHists->Add(fEventStat);
421 fHists->Add(fHistPt);
422 fHists->Add(fNContrNtracklets);
423 fHists->Add(fNContrNtracks);
424 fHists->Add(fCorruptedChunks);
425 fHists->Add(fCorruptedChunksAfter);
428 fHists->Add(fHistPtMC);
430 fHists->Add(fNmcNch);
431 fHists->Add(fPNmcNch);
432 fHists->Add(fNmcNchVtx);
433 fHists->Add(fNmcNchVtxStrangeCorr);
434 fHists->Add(fPNmcNchVtx);
436 fHists->Add(fNmcNchTracklet);
437 fHists->Add(fPNmcNchTracklet);
438 fHists->Add(fNmcNchVtxTracklet);
439 fHists->Add(fPNmcNchVtxTracklet);
441 fHists->Add(fChargedPi0);
442 fHists->Add(fVertexCheck);
443 fHists->Add(fPropagateDca);
445 for(Int_t i=0;i<8;i++){
446 fHists->Add(fMapSingleTrig[i]);
447 fHists->Add(fMapPair[i]);
448 fHists->Add(fMapEvent[i]);
449 fHists->Add(fMapAll[i]);
450 fHists->Add(fMapThree[i]);
451 fHists->Add(fVertexZ[i]);
453 fHists->Add(fNcharge[i]);
454 fHists->Add(fEta[i]);
455 fHists->Add(fPhi[i]);
456 fHists->Add(fDcaXY[i]);
457 fHists->Add(fDcaZ[i]);
458 fHists->Add(fPtSeed[i]);
459 fHists->Add(fEtaSeed[i]);
460 fHists->Add(fPhiSeed[i]);
461 fHists->Add(fPtOthers[i]);
462 fHists->Add(fEtaOthers[i]);
463 fHists->Add(fPhiOthers[i]);
464 fHists->Add(fPtEtaOthers[i]);
465 fHists->Add(fPhiEta[i]);
466 fHists->Add(fDPhiDEtaEventAxis[i]);
467 fHists->Add(fTriggerNch[i]);
468 fHists->Add(fTriggerNchSeeds[i]);
469 fHists->Add(fTriggerTracklet[i]);
470 fHists->Add(fNch07Nch[i]);
471 fHists->Add(fPNch07Nch[i]);
472 fHists->Add(fNch07Tracklet[i]);
473 fHists->Add(fNchTracklet[i]);
474 fHists->Add(fPNch07Tracklet[i]);
475 fHists->Add(fDPhiEventAxis[i]);
476 fHists->Add(fDPhi1DPhi2[i]);
483 //________________________________________________________________________
484 void AliAnalysisTaskMinijet::UserExec(Option_t *)
486 // Main function, called for each event
487 // Kinematics-only, ESD and AOD can be processed.
488 // Data is read (ReadEventESD, ReadEventAOD...) and then analysed (Analyse).
489 // - in case of MC with full detector simulation, all correction steps(0-5) can be processed
490 // - for Data, only step 5 is performed
491 // - for kinematics-only, only step 0 is processed
493 // - trigger - - vertex - - tracks - - multiplicity -
494 // step 7 = Triggered events, reconstructed accepted vertex, reconstructed tracks + strangness corr, reconstructed multiplicity+strangeCorr
495 // step 6 = Triggered events, reconstructed accepted vertex, reconstructed tracks with MC prop + stangess corr, reconstructed multiplicity+strangeCorr
496 // step 5 = Triggered events, reconstructed accepted vertex, reconstructed tracks, reconstructed multiplicity
497 // step 4 = Triggered events, reconstructed accepted vertex, reconstructed tracks with MC properties, reconstructed multiplicity
498 // step 3 = Triggered events, reconstructed accepted vertex, mc primary particles, reconstructed multiplicity
499 // step 2 = Triggered events, reconstructed accepted vertex, mc primary particles, true multiplicity
500 // step 1 = Triggered events, all mc primary particles, true multiplicity
501 // step 0 = All events, all mc primary particles, true multiplicity
504 if(fDebug) Printf("UserExec: Event starts");
506 AliVEvent *event = InputEvent();
508 Error("UserExec", "Could not retrieve event");
511 fBSign= event->GetMagneticField();
513 //get events, either ESD or AOD
514 fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
515 fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
520 vector<Short_t> charge;
521 vector<Float_t> strangenessWeight;
525 vector<Float_t> theta;
528 //number of accepted tracks and tracklets
529 Int_t ntracks = 0; //return value for reading functions for ESD and AOD
530 //Int_t ntracksRemove = 0; //return value for reading functions for ESD and AOD
531 vector<Int_t> nTracksTracklets; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc,
532 //for real nall=ncharged)
534 if(!fAODEvent && !fESDEvent)return;
537 Double_t centrality = 0;
538 if (fCentralityMethod.Length() > 0)
540 AliCentrality *centralityObj = 0;
542 centralityObj = fAODEvent->GetHeader()->GetCentralityP();
544 centralityObj = fESDEvent->GetCentrality();
547 centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
548 AliInfo(Form("Centrality is %f", centrality));
552 Printf("WARNING: Centrality object is 0");
557 // SDD check for LHC11a
563 const AliMultiplicity *mul = fESDEvent->GetMultiplicity();
564 Int_t nClu3 = mul->GetNumberOfITSClusters(2);
565 Int_t nClu4 = mul->GetNumberOfITSClusters(3);
566 if(nClu3==0 && nClu4==0) return;
568 else if (fSelOption==1){
569 TString trcl = fESDEvent->GetFiredTriggerClasses().Data();
570 if (!(trcl.Contains("CINT1-B-NOPF-ALLNOTRD"))) return;
577 Bool_t useEvent = false;
578 Int_t nTracks = fAODEvent->GetNTracks();
579 for(Int_t itrack=0; itrack<nTracks; itrack++) {
580 AliAODTrack * track = fAODEvent->GetTrack(itrack);
581 if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
586 if (!useEvent) return;
588 else if(fSelOption==1){
589 TString trcl = fAODEvent->GetFiredTriggerClasses().Data();
590 if (!(trcl.Contains("CINT1-B-NOPF-ALLNOTRD"))) return;
596 fNMcPrimAccept=0;// number of accepted primaries
597 fNRecAccept=0; // number of accepted tracks
598 fNRecAcceptStrangeCorr=0; // number of accepted tracks + strangeness correction
599 fNMcPrimAcceptTracklet=0;// number of accepted primaries (no pt cut)
600 fNRecAcceptTracklet=0; // number of accepted tracklets
602 // instead of task->SelectCollisionCandidate(mask) in AddTask macro
603 Bool_t isSelected = (((AliInputEventHandler*)
604 (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
605 ->IsEventSelected() & fTriggerType);
609 Printf("IsSelected = %d", isSelected);
610 Printf("CheckEvent(true)= %d", CheckEvent(true));
611 Printf("CheckEvent(false)= %d", CheckEvent(false));
614 fEventStat->Fill(0);//all events
617 if(isSelected){ // has offline trigger
619 fEventStat->Fill(1);//triggered event
621 if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec
623 fEventStat->Fill(2);//triggered event with vertex
626 //step 5 = TrigVtxRecNrec
629 if(fESDEvent) ntracks = ReadEventESD(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 5);
630 else if(fAODEvent)ntracks = ReadEventAOD(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 5);
631 else AliInfo("Fatal Error");
633 if (fCentralityMethod.Length() > 0)
634 ntracks = TMath::Nint(centrality);
637 if(pt.size()){ //(internally ntracks=fNRecAccept)
638 fEventStat->Fill(3);//triggered event with vertex and one reconstructed track in acceptance
639 Analyse(pt, eta, phi, charge, strangenessWeight, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//step 5 = TrigVtxRecNrec
642 if(fCorrStrangeness){
643 //step 7 = TrigVtxRecNrecStrangeCorr
646 if(fESDEvent) ntracks = ReadEventESD(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 7);//stagness version not yet implemented
647 else if(fAODEvent)ntracks = ReadEventAOD(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 7);
648 else AliInfo("Fatal Error");
650 if (fCentralityMethod.Length() > 0)
651 ntracks = TMath::Nint(centrality);
654 if(pt.size()){ //(internally ntracks=fNRecAccept)
655 Analyse(pt, eta, phi, charge, strangenessWeight, TMath::Nint(fNRecAcceptStrangeCorr), nTracksTracklets[1], nTracksTracklets[2], 7);//step 7 = TrigVtxRecNrecStrangeCorr
660 // step 4 = TrigVtxRecMcPropNrec
663 if(fESDEvent) ntracks = ReadEventESDRecMcProp(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 4);
664 else if(fAODEvent) ntracks = ReadEventAODRecMcProp(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 4);
665 else AliInfo("Fatal Error");
667 if (fCentralityMethod.Length() > 0)
668 ntracks = TMath::Nint(centrality);
671 if(pt.size()){//(internally ntracks=fNRecAccept)
672 Analyse(pt, eta, phi, charge,strangenessWeight, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4); //step 4 = TrigVtxRecMcPropNrec
676 if(fCorrStrangeness){
677 // step 6 = TrigVtxRecMcPropNrecStrangeCorr
680 if(fESDEvent) ntracks = ReadEventESDRecMcProp(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 6);//stagness version not yet implemented
681 else if(fAODEvent) ntracks = ReadEventAODRecMcProp(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 6);
682 else AliInfo("Fatal Error");
684 if (fCentralityMethod.Length() > 0)
685 ntracks = TMath::Nint(centrality);
688 if(pt.size()){//(internally ntracks=fNRecAccept)
689 Analyse(pt, eta, phi, charge, strangenessWeight, TMath::Nint(fNRecAcceptStrangeCorr), nTracksTracklets[1], nTracksTracklets[2], 6); //step 6 = TrigVtxRecMcPropNrecStrangeCorr
692 // step 3 = TrigVtxMcNrec
695 if(fESDEvent) ntracks = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 3);
696 else if(fAODEvent) ntracks = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 3);
697 else AliInfo("Fatal Error");
699 if (fCentralityMethod.Length() > 0){
700 fNRecAccept = TMath::Nint(centrality);
701 fNMcPrimAccept = TMath::Nint(centrality);
706 Analyse(pt, eta, phi, charge, strangenessWeight, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 3); //step 3 = TrigVtxMcNrec
707 Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 2); //step 2 = TrigVtxMcNmc
713 }//check event (true)
716 if(fUseMC && !fMcOnly){
718 fNMcPrimAccept=0;// number of accepted primaries
719 fNRecAccept=0; // number of accepted tracks
720 fNMcPrimAcceptTracklet=0;// number of accepted primaries (no pt cut)
721 fNRecAcceptTracklet=0; // number of accepted tracklets
723 if(CheckEvent(false)){// all events, with and without reconstucted vertex
724 if(fESDEvent) ntracks = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 1);//read tracks
725 else if(fAODEvent) ntracks = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 1);//read tracks
726 else AliInfo("Fatal Error");
728 if (fCentralityMethod.Length() > 0)
729 fNMcPrimAccept = TMath::Nint(centrality);
734 Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1); // step 1 = TrigAllMcNmc
736 Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0); //first part of step 0 // step 0 = AllAllMcNmc
744 else { // not selected by physics selection task = not triggered
745 if(fUseMC && !fMcOnly){
746 if(CheckEvent(false)){
749 if(fESDEvent) ntracks = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
750 else if(fAODEvent) ntracks = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
751 else AliInfo("Fatal Error");
753 if (fCentralityMethod.Length() > 0)
754 fNMcPrimAccept = TMath::Nint(centrality);
758 Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0); //second part of step 0 // step 0 = AllAllMcNmc
766 if(fMode==0) ntracks = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
767 else if (fMode==1) ntracks = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
769 if (fCentralityMethod.Length() > 0)
770 fNMcPrimAccept = TMath::Nint(centrality);
774 Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);
780 //________________________________________________________________________
781 Int_t AliAnalysisTaskMinijet::ReadEventESD( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
782 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
783 vector<Float_t> &strangeArray,
784 vector<Int_t> &nTracksTracklets, const Int_t step)
786 // gives back the number of esd tracks and pointer to arrays with track
787 // properties (pt, eta, phi)
788 // Uses TPC tracks with SPD vertex from now on
794 strangeArray.clear();
795 nTracksTracklets.clear();
797 const AliESDVertex* vtxSPD = fESDEvent->GetPrimaryVertexSPD(); // uses track or SPD vertexer
798 fVertexZ[step]->Fill(vtxSPD->GetZ());
800 // Retreive the number of all tracks for this event.
801 Int_t ntracks = fESDEvent->GetNumberOfTracks();
802 if(fDebug>1) Printf("all ESD tracks: %d", ntracks);
804 //first loop to check how many tracks are accepted
806 Int_t nAcceptedTracks=0;
807 //Float_t nAcceptedTracksStrange=0;
808 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
810 AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
812 Error("ReadEventESD", "Could not receive track %d", iTracks);
816 // use TPC only tracks with non default SPD vertex
817 AliESDtrack *track = AliESDtrackCuts::
818 GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
820 if(!fCuts->AcceptTrack(track)) {
822 continue;// apply TPC track cuts before constrain to SPD vertex
825 // only constrain tracks above threshold
826 AliExternalTrackParam exParam;
827 // take the B-field from the ESD, no 3D fieldMap available at this point
828 Bool_t relate = false;
829 relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
835 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
836 exParam.GetCovariance());
840 continue;// only if tracks have pt<=0
843 if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax) {
844 ptArray.push_back(track->Pt());
845 etaArray.push_back(track->Eta());
846 phiArray.push_back(track->Phi());
847 chargeArray.push_back(track->Charge());
848 strangeArray.push_back(1);
850 fHistPt->Fill(track->Pt());
853 // TPC only track memory needs to be cleaned up
860 if(nAcceptedTracks==0) return 0;
863 Int_t ntrackletsAccept=0;
864 AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
865 Int_t ntracklets = mult->GetNumberOfTracklets();
866 for(Int_t i=0;i< ntracklets;i++){
867 if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
871 nTracksTracklets.push_back(nAcceptedTracks);
872 nTracksTracklets.push_back(ntrackletsAccept);
873 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
874 //where here also neutral particles are counted.
877 fVzEvent=vtxSPD->GetZ(); // needed for correction map
878 if(step==5 || step ==2) {
879 fNRecAccept=nAcceptedTracks;
880 fNRecAcceptTracklet=ntrackletsAccept;
887 //________________________________________________________________________
888 Int_t AliAnalysisTaskMinijet::ReadEventESDRecMcProp( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
889 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
890 vector<Float_t> &strangeArray,
891 vector<Int_t> &nTracksTracklets, const Int_t step)
893 // gives back the number of esd tracks and pointer to arrays with track
894 // properties (pt, eta, phi) of mc particles if available
895 // Uses TPC tracks with SPD vertex from now on
901 strangeArray.clear();
902 nTracksTracklets.clear();
905 AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
907 Error("ReadEventESDRecMcProp", "Could not retrieve MC event");
910 AliStack* stack = MCEvent()->Stack();
914 // Retreive the number of all tracks for this event.
915 Int_t ntracks = fESDEvent->GetNumberOfTracks();
916 if(fDebug>1)Printf("all ESD tracks: %d", ntracks);
918 const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
919 fVertexZ[step]->Fill(vtxSPD->GetZ());
922 Int_t nAcceptedTracks=0;
923 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
925 AliVParticle *vtrack = fESDEvent->GetTrack(iTracks);
926 AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
928 Error("ReadEventESDRecMcProp", "Could not receive track %d", iTracks);
932 // use TPC only tracks with non default SPD vertex
933 AliESDtrack *track = AliESDtrackCuts::
934 GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
936 if(!fCuts->AcceptTrack(track)) {
938 continue;// apply TPC track cuts before constrain to SPD vertex
941 // only constrain tracks above threshold
942 AliExternalTrackParam exParam;
943 // take the B-field from the ESD, no 3D fieldMap available at this point
944 Bool_t relate = false;
945 relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
951 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
952 exParam.GetCovariance());
956 continue;// only if tracks have pt<=0
959 //count tracks, if available, use mc particle properties
960 if(vtrack->GetLabel()<=0){
961 if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax){
962 ptArray.push_back(track->Pt());
963 etaArray.push_back(track->Eta());
964 phiArray.push_back(track->Phi());
965 chargeArray.push_back(track->Charge());
966 strangeArray.push_back(1);
971 TParticle *partOfTrack = stack->Particle(vtrack->GetLabel());
972 if (TMath::Abs(partOfTrack->Eta())<fEtaCut && partOfTrack->Pt()>fPtMin && partOfTrack->Pt()<fPtMax) {
973 ptArray.push_back(partOfTrack->Pt());
974 etaArray.push_back(partOfTrack->Eta());
975 phiArray.push_back(partOfTrack->Phi());
976 chargeArray.push_back(vtrack->Charge());
977 strangeArray.push_back(1);
982 // TPC only track memory needs to be cleaned up
988 if(nAcceptedTracks==0) return 0;
991 Int_t ntrackletsAccept=0;
992 AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
993 Int_t ntracklets = mult->GetNumberOfTracklets();
994 for(Int_t i=0;i< ntracklets;i++){
995 if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
1000 nTracksTracklets.push_back(nAcceptedTracks);
1001 nTracksTracklets.push_back(ntrackletsAccept);
1002 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1003 //where here also neutral particles are counted.
1006 //get mc vertex for correction maps
1007 AliGenEventHeader* header = MCEvent()->GenEventHeader();
1009 header->PrimaryVertex(mcV);
1012 return fNRecAccept; // give back reconstructed value
1020 //________________________________________________________________________
1021 Int_t AliAnalysisTaskMinijet::ReadEventESDMC(vector<Float_t> &ptArray, vector<Float_t> &etaArray,
1022 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1023 vector<Float_t> &strangeArray,
1024 vector<Int_t> &nTracksTracklets, const Int_t step)
1026 // gives back the number of charged prim MC particle and pointer to arrays
1027 // with particle properties (pt, eta, phi)
1032 chargeArray.clear();
1033 strangeArray.clear();
1034 nTracksTracklets.clear();
1038 AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
1040 Error("ReadEventESDMC", "Could not retrieve MC event");
1044 AliStack* stack = MCEvent()->Stack();
1045 if(!stack) return 0;
1047 Int_t ntracks = mcEvent->GetNumberOfTracks();
1048 if(fDebug>1)Printf("MC particles: %d", ntracks);
1051 AliGenEventHeader* header = MCEvent()->GenEventHeader();
1055 header->PrimaryVertex(mcV);
1058 fVertexZ[0]->Fill(vzMC);
1060 fVertexZ[step]->Fill(vzMC);
1063 //----------------------------------
1064 //first track loop to check how many chared primary tracks are available
1065 Int_t nChargedPrimaries=0;
1066 Int_t nAllPrimaries=0;
1068 Int_t nPseudoTracklets=0;
1069 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1070 AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
1072 Error("ReadEventESDMC", "Could not receive track %d", iTracks);
1077 if(//count also charged particles in case of fSelectParticles==2 (only neutral)
1078 !SelectParticlePlusCharged(
1081 stack->IsPhysicalPrimary(track->Label())
1086 //count number of pseudo tracklets
1087 if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.0) ++nPseudoTracklets; //0.035
1088 //same cuts as on ESDtracks
1089 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1091 //count all primaries
1093 //count charged primaries
1094 if (track->Charge() != 0) ++nChargedPrimaries;
1096 if(fDebug>2) Printf("PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
1100 //----------------------------------
1103 Printf("All in acceptance=%d",nAllPrimaries);
1104 Printf("Charged in acceptance =%d",nChargedPrimaries);
1107 fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
1109 if(nAllPrimaries==0) return 0;
1110 if(nChargedPrimaries==0) return 0;
1113 //Int_t iChargedPiK=0;
1114 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1115 AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
1117 Error("ReadEventESDMC", "Could not receive track %d", iTracks);
1124 stack->IsPhysicalPrimary(track->Label())
1130 //same cuts as on ESDtracks
1131 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1133 if(fDebug>2) Printf("After: PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
1136 fHistPtMC->Fill(track->Pt());
1137 //fills arrays with track properties
1138 ptArray.push_back(track->Pt());
1139 etaArray.push_back(track->Eta());
1140 phiArray.push_back(track->Phi());
1141 chargeArray.push_back(track->Charge());
1142 strangeArray.push_back(1);
1146 nTracksTracklets.push_back(nChargedPrimaries);
1147 nTracksTracklets.push_back(nPseudoTracklets);
1148 if(fSelectParticles!=2){
1149 nTracksTracklets.push_back(nAllPrimaries);
1152 nTracksTracklets.push_back(nAllPrimaries-nChargedPrimaries); // only neutral
1155 fNMcPrimAccept = nChargedPrimaries;
1156 fNMcPrimAcceptTracklet = nPseudoTracklets;
1159 fNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1160 fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1161 fNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1162 fPNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1165 fNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1166 fPNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1167 fNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1168 fPNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1176 //________________________________________________________________________
1177 Int_t AliAnalysisTaskMinijet::ReadEventAOD( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
1178 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1179 vector<Float_t> &strangeArray,
1180 vector<Int_t> &nTracksTracklets, const Int_t step)
1182 // gives back the number of AOD tracks and pointer to arrays with track
1183 // properties (pt, eta, phi)
1188 chargeArray.clear();
1189 strangeArray.clear();
1190 nTracksTracklets.clear();
1192 TClonesArray *mcArray=0x0;
1193 if(fAnalysePrimOnly || (fCorrStrangeness && fUseMC)){
1194 mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
1198 AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1199 Double_t vzAOD=vertex->GetZ();
1200 fVertexZ[step]->Fill(vzAOD);
1202 // Retreive the number of tracks for this event.
1203 Int_t ntracks = fAODEvent->GetNumberOfTracks();
1204 if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1207 Int_t nAcceptedTracks=0;
1208 Float_t nAcceptedTracksStrange=0;
1209 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1210 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
1212 Error("ReadEventAOD", "Could not receive track %d", iTracks);
1216 AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1218 //use only tracks from primaries
1219 if(fAnalysePrimOnly){
1220 if(vtrack->GetLabel()<=0)continue;
1221 if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1224 Double_t save= track->Pt();
1225 Double_t d0rphiz[2],covd0[3];
1226 Bool_t isDca= track->PropagateToDCA(fAODEvent->GetPrimaryVertex(),fAODEvent->GetMagneticField(),9999.,d0rphiz,covd0);
1227 fPropagateDca->Fill(Int_t(isDca));
1228 if(TMath::Abs(save - track->Pt())>1e-6) Printf("Before pt=%f, After pt=%f",save, track->Pt());
1230 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
1231 && track->Pt()>fPtMin && track->Pt()<fPtMax){
1235 // Printf("dca= %f", track->DCA());
1236 //save track properties in vector
1237 ptArray.push_back(track->Pt());
1238 etaArray.push_back(track->Eta());
1239 phiArray.push_back(track->Phi());
1240 chargeArray.push_back(track->Charge());
1241 fHistPt->Fill(track->Pt());
1244 //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1246 if(fUseMC && fCorrStrangeness && step==7){
1247 if(vtrack->GetLabel()>0){
1248 AliAODMCParticle* mcprong =(AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1250 Int_t labmom = mcprong->GetMother();
1252 AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1255 pdgMother = mcmother->GetPdgCode();
1256 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1257 if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1258 else factor=1./0.6;// values from strangeness publication
1260 if(TMath::Abs(pdgMother)==3122) { //Lambda
1261 factor=1./0.25; // values from strangeness publication
1268 nAcceptedTracksStrange+=factor;
1269 strangeArray.push_back(factor);
1270 fDcaXY[step]->Fill(d0rphiz[0], factor);
1271 fDcaZ[step]->Fill(d0rphiz[0], factor);
1275 //need to check this option for MC
1276 if(nAcceptedTracks==0) return 0;
1280 Int_t ntrackletsAccept=0;
1281 AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1282 for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1283 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut){
1289 nTracksTracklets.push_back(nAcceptedTracks);
1290 nTracksTracklets.push_back(ntrackletsAccept);
1291 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1292 //where here also neutral particles are counted.
1296 if(step==5 || step==2){
1297 fNRecAccept = nAcceptedTracks; // needed for MC case //step5 = TrigVtxRecNrec
1298 fNRecAcceptTracklet = ntrackletsAccept; // needed for MC case //step5 = TrigVtxRecNrec
1300 if(step==7)fNRecAcceptStrangeCorr = nAcceptedTracksStrange;
1302 return fNRecAccept; // at the moment, always return reconstructed multiplicity
1306 //________________________________________________________________________
1307 Int_t AliAnalysisTaskMinijet::ReadEventAODRecMcProp( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
1308 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1309 vector<Float_t> &strangeArray,
1310 vector<Int_t> &nTracksTracklets, const Int_t step)
1312 // gives back the number of AOD tracks and pointer to arrays with track
1313 // properties (pt, eta, phi)
1319 chargeArray.clear();
1320 strangeArray.clear();
1321 nTracksTracklets.clear();
1324 // Retreive the number of tracks for this event.
1325 Int_t ntracks = fAODEvent->GetNumberOfTracks();
1326 if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1329 //get array of mc particles
1330 TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1331 FindListObject(AliAODMCParticle::StdBranchName());
1333 AliInfo("No MC particle branch found");
1337 AliAODVertex* vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1338 Double_t vzAOD=vtx->GetZ();
1339 fVertexZ[step]->Fill(vzAOD);
1341 Int_t nAcceptedTracks=0;
1342 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1343 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
1345 AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1348 Error("ReadEventAODRecMcProp", "Could not receive track %d", iTracks);
1352 //use only tracks from primaries
1353 if(fAnalysePrimOnly){
1354 if(vtrack->GetLabel()<=0)continue;
1355 if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1358 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut &&
1359 track->Pt()>fPtMin && track->Pt()<fPtMax){
1364 //save track properties in vector
1365 if(vtrack->GetLabel()<=0){ //fake tracks before "label<0", but crash in AOD079 // what is the meaning of label 0
1366 // Printf("Fake track");
1368 ptArray.push_back(track->Pt());
1369 etaArray.push_back(track->Eta());
1370 phiArray.push_back(track->Phi());
1371 chargeArray.push_back(track->Charge());
1374 else{//mc properties
1375 AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1376 if(!partOfTrack) Printf("label=%d", vtrack->GetLabel());
1378 ptArray.push_back(partOfTrack->Pt());
1379 etaArray.push_back(partOfTrack->Eta());
1380 phiArray.push_back(partOfTrack->Phi());
1381 chargeArray.push_back(vtrack->Charge());//partOfTrack?
1383 //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1384 if(fUseMC && fCorrStrangeness && step==6){
1385 Int_t labmom = partOfTrack->GetMother();
1387 AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1390 pdgMother = mcmother->GetPdgCode();
1391 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1392 if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1393 else factor=1./0.6;// values from strangeness publication
1395 if(TMath::Abs(pdgMother)==3122) { //Lambda
1396 factor=1./0.25; // values from strangeness publication
1403 strangeArray.push_back(factor);
1407 //need to check this option for MC
1408 if(nAcceptedTracks==0) return 0;
1411 Int_t ntrackletsAccept=0;
1412 AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1413 for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1414 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut ){
1420 nTracksTracklets.push_back(nAcceptedTracks);
1421 nTracksTracklets.push_back(ntrackletsAccept);
1422 nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1423 //where here also neutral particles are counted.
1427 AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1428 FindListObject(AliAODMCHeader::StdBranchName());
1429 Float_t vzMC = aodMCheader->GetVtxZ();
1432 return fNRecAccept;//this is the rec value from step 5
1438 //________________________________________________________________________
1439 Int_t AliAnalysisTaskMinijet::ReadEventAODMC( vector<Float_t> &ptArray, vector<Float_t> &etaArray,
1440 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1441 vector<Float_t> &strangeArray,
1442 vector<Int_t> &nTracksTracklets, const Int_t step)
1444 // gives back the number of AOD MC particles and pointer to arrays with particle
1445 // properties (pt, eta, phi)
1450 chargeArray.clear();
1451 strangeArray.clear();
1452 nTracksTracklets.clear();
1455 AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1456 FindListObject(AliAODMCHeader::StdBranchName());
1457 Float_t vzMC = aodMCheader->GetVtxZ();
1459 fVertexZ[0]->Fill(vzMC);
1461 fVertexZ[step]->Fill(vzMC);
1463 //retreive MC particles from event
1464 TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1465 FindListObject(AliAODMCParticle::StdBranchName());
1467 AliInfo("No MC particle branch found");
1471 Int_t ntracks = mcArray->GetEntriesFast();
1472 if(fDebug>1)Printf("MC particles: %d", ntracks);
1475 // Track loop: chek how many particles will be accepted
1477 Int_t nChargedPrim=0;
1479 Int_t nPseudoTracklets=0;
1480 for (Int_t it = 0; it < ntracks; it++) {
1481 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1483 Error("ReadEventAODMC", "Could not receive particle %d", it);
1487 if(!SelectParticlePlusCharged(
1490 track->IsPhysicalPrimary()
1495 if(TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.0)++nPseudoTracklets; //0.035
1496 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1499 if(track->Charge()!=0) nChargedPrim++;
1504 if(nAllPrim==0) return 0;
1505 if(nChargedPrim==0) return 0;
1507 //generate array with size of number of accepted tracks
1508 fChargedPi0->Fill(nAllPrim,nChargedPrim);
1511 // Track loop: fill arrays for accepted tracks
1512 // Int_t iChargedPrim=0;
1513 for (Int_t it = 0; it < ntracks; it++) {
1514 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1516 Error("ReadEventAODMC", "Could not receive particle %d", it);
1523 track->IsPhysicalPrimary()
1528 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1530 fHistPtMC->Fill(track->Pt());
1531 ptArray.push_back(track->Pt());
1532 etaArray.push_back(track->Eta());
1533 phiArray.push_back(track->Phi());
1534 chargeArray.push_back(track->Charge());
1535 strangeArray.push_back(1);
1538 nTracksTracklets.push_back(nChargedPrim);
1539 nTracksTracklets.push_back(nPseudoTracklets);
1540 if(fSelectParticles!=2){
1541 nTracksTracklets.push_back(nAllPrim);
1544 nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral
1550 fNMcPrimAccept = nChargedPrim;
1551 fNMcPrimAcceptTracklet = nPseudoTracklets;
1553 if(step==1){ // step 1 = Trig All Mc Nmc
1554 fNmcNch->Fill( fNMcPrimAccept,fNRecAccept);
1555 fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1556 fNmcNchTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1557 fPNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1559 if(step==3){ // step 3 = Trig vtx Mc
1560 fNmcNchVtx->Fill( fNMcPrimAccept,fNRecAccept);
1561 fNmcNchVtxStrangeCorr->Fill( fNMcPrimAccept,fNRecAcceptStrangeCorr);
1562 fPNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1563 fNmcNchVtxTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1564 fPNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1566 return fNRecAccept; // rec value from step 5 or step 2
1571 //________________________________________________________________________
1572 void AliAnalysisTaskMinijet::Analyse(const vector<Float_t> &pt,
1573 const vector<Float_t> &eta,
1574 const vector<Float_t> &phi,
1575 const vector<Short_t> &charge,
1576 const vector<Float_t> &strangeWeight,
1577 const Int_t ntracksCharged,
1578 const Int_t ntracklets,
1583 // analyse track properties (comming from either ESDs or AODs) in order to compute
1584 // mini jet activity (chared tracks) as function of charged multiplicity
1589 Printf("Analysis Step=%d", step);
1591 Printf("nAll=%d",nAll);
1592 Printf("nCharged=%d",ntracksCharged);
1593 for (Int_t i = 0; i < nAll; i++) {
1594 Printf("pt[%d]=%f",i,pt[i]);
1599 //calculation of mean pt for all tracks in case of step==0
1600 if(step==5 || step ==2){ // rec step
1602 Double_t leadingPt=0.;
1603 for (Int_t i = 0; i < nAll; i++) {
1604 if(pt[i]<0.01)continue;
1606 if(leadingPt<pt[i])leadingPt=pt[i];
1610 fLeadingPtRec=leadingPt;
1613 Double_t propEvent[4] = {ntracksCharged,fVzEvent,fMeanPtRec, fLeadingPtRec}; //vz: {rec, mc, mc}, meanPt and Nrec is always rec value
1614 fMapEvent[step]->Fill(propEvent);
1616 fNcharge[step]->Fill(ntracksCharged);
1618 Float_t ptEventAxis=0; // pt event axis
1619 Float_t etaEventAxis=0; // eta event axis
1620 Float_t phiEventAxis=0; // phi event axis
1621 Short_t chargeEventAxis=0; // charge event axis
1622 Float_t strangeWeightEventAxis=0; // strange weight event axis
1624 Float_t ptOthers = 0; // pt others // for all other tracks around event axis -> see loop
1625 Float_t etaOthers = 0; // eta others
1626 Float_t phiOthers = 0; // phi others
1627 Short_t chargeOthers = 0; // charge others
1628 Float_t strangeWeightOthers = 0; // strange weight others
1630 Int_t *pindexInnerEta = new Int_t [nAll+1];
1631 Float_t *ptInnerEta = new Float_t[nAll+1];
1635 for (Int_t i = 0; i < nAll; i++) {
1637 if(pt[i]<0.01)continue;
1639 //fill single particle correction for first step of pair correction
1640 Double_t propAll[3] = {eta[i],pt[i],ntracksCharged };
1641 fMapAll[step]->Fill(propAll, strangeWeight[i]);
1644 //filling of simple check plots
1645 if(pt[i]<0.7) continue;
1646 fPt[step] -> Fill( pt[i]);
1647 fEta[step] -> Fill(eta[i]);
1648 fPhi[step] -> Fill(phi[i]);
1649 fPhiEta[step]-> Fill(phi[i], eta[i]);
1651 pindexInnerEta[i]=0; //set all values to zero
1652 //fill new array for eta check
1653 ptInnerEta[i]= pt[i];
1659 // define event axis: leading or random track (with pt>fTriggerPtCut)
1660 // ---------------------------------------
1661 Int_t highPtTracks=0;
1662 Int_t highPtTracksInnerEta=0;
1663 // Double_t highPtTracksInnerEtaCorr=0;
1666 //count high pt tracks and high pt tracks in acceptance for seeds
1667 for(Int_t i=0;i<nAll;i++){
1669 if(pt[i]<0.01)continue;
1671 if(TMath::Abs(eta[i])<0.9){
1675 if(pt[i]>fTriggerPtCut) {
1679 // seed should be place in middle of acceptance, that complete cone is in acceptance
1680 if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
1682 // Printf("eta=%f", eta[i]);
1683 highPtTracksInnerEta++;
1692 //sort array in order to get highest pt tracks first
1693 //index can be computed with pindexInnerEta[number]
1694 if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1696 // plot of multiplicity distributions
1697 fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1698 fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1701 fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1702 fNchTracklet[step]->Fill(ntracklets, ntracksCharged);
1703 fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1706 //analysis can only be performed with event axis, defined by high pt track
1709 if(highPtTracks>0 && highPtTracksInnerEta>0){
1711 // build pairs in two track loops
1712 // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
1713 for(Int_t axis=0;(axis<nAll) && (pt[pindexInnerEta[axis]]>=fTriggerPtCut); axis++){
1715 //EventAxisRandom track properties
1716 ptEventAxis = pt [pindexInnerEta[axis]];
1717 etaEventAxis = eta[pindexInnerEta[axis]];
1718 phiEventAxis = phi[pindexInnerEta[axis]];
1719 chargeEventAxis = charge[pindexInnerEta[axis]];
1720 strangeWeightEventAxis = strangeWeight[pindexInnerEta[axis]];
1721 fPtSeed[step] -> Fill( ptEventAxis);
1722 fEtaSeed[step] -> Fill(etaEventAxis);
1723 fPhiSeed[step] -> Fill(phiEventAxis);
1726 Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged };
1727 fMapSingleTrig[step]->Fill(prop, strangeWeightEventAxis);
1730 for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) {
1732 if(pt[pindexInnerEta[iTrack]]<fAssociatePtCut) continue;
1734 if(fSelectParticlesAssoc==1){
1735 if(charge[pindexInnerEta[iTrack]]==0)continue;
1737 if(fSelectParticlesAssoc==2){
1738 if(charge[pindexInnerEta[iTrack]]!=0)continue;
1742 ptOthers = pt [pindexInnerEta[iTrack]];
1743 etaOthers = eta[pindexInnerEta[iTrack]];
1744 phiOthers = phi[pindexInnerEta[iTrack]];
1745 chargeOthers = charge[pindexInnerEta[iTrack]];
1746 strangeWeightOthers = strangeWeight[pindexInnerEta[iTrack]];
1749 //plot only properties of tracks with pt>ptassoc
1750 fPtOthers[step] -> Fill( ptOthers);
1751 fEtaOthers[step] -> Fill(etaOthers);
1752 fPhiOthers[step] -> Fill(phiOthers);
1753 fPtEtaOthers[step] -> Fill(ptOthers, etaOthers);
1755 // if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]);
1757 Float_t dPhi = phiOthers-phiEventAxis;
1758 if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1759 else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1760 Float_t dEta=etaOthers-etaEventAxis;
1763 fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta, strangeWeightEventAxis*strangeWeightOthers);
1764 fDPhiEventAxis[step]->Fill(dPhi, strangeWeightEventAxis*strangeWeightOthers);
1767 if(ptEventAxis< 0.4 || ptEventAxis > 100) AliInfo("particles out of range pt");
1768 if(ntracksCharged<-1 || ntracksCharged>1500) AliInfo("particles out of range ncharge");
1769 if(TMath::Abs(dEta)>2*fEtaCut) {
1770 AliInfo("particles out of range dEta");
1771 AliInfo(Form("eta1=%f, eta2=%f", etaOthers, etaEventAxis));
1772 AliInfo(Form("step=%d",step));
1774 if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){
1775 AliInfo(Form("particles out of range dPhi"));
1776 AliInfo(Form("phi1=%f, phi2=%f", phiOthers, phiEventAxis));
1779 Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers);
1781 Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, isLikeSign };
1782 fMapPair[step]->Fill(prop6, strangeWeightEventAxis*strangeWeightOthers);
1784 //thrid track loop (Andreas: three particle auto-correlations)
1785 if(fThreeParticleCorr){
1786 for (Int_t third = iTrack+1; third < nAll; third++) {
1787 if(pt[pindexInnerEta[iTrack]]<fTriggerPtCut) continue;
1788 if(pt[pindexInnerEta[third]]<fTriggerPtCut) continue;
1791 Float_t dPhi1 = phiEventAxis - phiOthers;
1792 if(dPhi1>1.5*TMath::Pi()) dPhi1 = dPhi1-2*TMath::Pi();
1793 else if(dPhi1<-0.5*TMath::Pi())dPhi1=dPhi1+2*TMath::Pi();
1795 Float_t phiThird = phi[pindexInnerEta[third]];
1796 Float_t strangeWeightThird = strangeWeight[pindexInnerEta[third]];
1799 Float_t dPhi2 = phiEventAxis - phiThird;
1800 if(dPhi2>1.5*TMath::Pi()) dPhi2 = dPhi2-2*TMath::Pi();
1801 else if(dPhi2<-0.5*TMath::Pi())dPhi2=dPhi2+2*TMath::Pi();
1803 fDPhi1DPhi2[step]-> Fill(dPhi1, dPhi2, strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird);
1804 Double_t propThree[3] = {dPhi1,dPhi2,ntracksCharged};
1805 fMapThree[step]->Fill(propThree,strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird );
1808 }// end of three particle correlation loop
1810 }// if fThreeParticleCorr is set to true
1812 }// end of inner track loop
1814 } //end of outer track loop
1816 fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1817 fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1818 fTriggerTracklet[step]->Fill(ntracklets);
1821 }//if there is at least one high pt track
1824 if(pindexInnerEta){// clean up array memory used for TMath::Sort
1825 delete[] pindexInnerEta;
1829 if(ptInnerEta){// clean up array memory used for TMath::Sort
1830 delete[] ptInnerEta;
1838 //________________________________________________________________________
1839 void AliAnalysisTaskMinijet::Terminate(Option_t*)
1841 //terminate function is called at the end
1842 //can be used to draw histograms etc.
1847 //________________________________________________________________________
1848 Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, Bool_t prim)
1850 //selection of mc particle
1851 //fSelectParticles=0: use charged primaries and pi0 and k0
1852 //fSelectParticles=1: use only charged primaries
1853 //fSelectParticles=2: use only pi0 and k0
1855 if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1857 if(!(pdg==111||pdg==130||pdg==310))
1866 else if(fSelectParticles==1){
1867 if (charge==0 || !prim){
1873 AliInfo("Error: wrong selection of charged/pi0/k0");
1881 //________________________________________________________________________
1882 Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
1884 //selection of mc particle
1885 //fSelectParticles=0: use charged primaries and pi0 and k0
1886 //fSelectParticles=1: use only charged primaries
1887 //fSelectParticles=2: use only pi0 and k0
1889 if(fSelectParticles==0){
1891 if(!(pdg==111||pdg==130||pdg==310))
1899 else if (fSelectParticles==1){
1900 if (charge==0 || !prim){
1904 else if(fSelectParticles==2){
1905 if(!(pdg==111||pdg==130||pdg==310))
1913 //________________________________________________________________________
1914 Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex)
1916 // This function tests the quality of an event (ESD/AOD) (rec and/or mc part)
1917 // recVertex=false: check if Mc events and stack is there, Nmc>0
1918 // recVertex=false: " + check if there is a good, reconstructed SPD vertex
1919 // defined by |z|<fVertexCut(10cm), Contributer>0, no PileUpFromSPD(3,0,8)
1927 AliMCEvent *mcEvente = (AliMCEvent*) MCEvent();
1929 Error("CheckEvent", "Could not retrieve MC event");
1934 AliStack* stackg = MCEvent()->Stack();
1935 if(!stackg) return false;
1936 Int_t ntracksg = mcEvente->GetNumberOfTracks();
1937 if(ntracksg<0) return false;
1940 AliGenEventHeader* headerg = MCEvent()->GenEventHeader();
1942 headerg->PrimaryVertex(mcVg);
1943 //if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 &&
1944 // TMath::Abs(mcVg[0])<1e-8) return false;
1945 Float_t vzMCg = mcVg[2];
1946 if(TMath::Abs(vzMCg)>fVertexZCut) return false;
1951 if(recVertex==true){
1952 if(fESDEvent->IsPileupFromSPD(3,0,8))return false;
1955 const AliESDVertex* vertexESDg = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
1956 if(!vertexESDg) return false;
1957 fVertexCheck->Fill(vertexESDg->GetNContributors());
1958 if(vertexESDg->GetNContributors()<=0)return false;
1959 Float_t fVzg= vertexESDg->GetZ();
1960 if(TMath::Abs(fVzg)>fVertexZCut) return false;
1963 const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
1964 if(!vtxSPD) return false;
1965 if(vtxSPD->GetNContributors()<=0)return false;
1966 Float_t fVzSPD= vtxSPD->GetZ();
1967 if(TMath::Abs(fVzSPD)>fVertexZCut) return false;
1974 else if(fMode==1){ //aod
1978 //retreive MC particles from event
1979 TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1980 FindListObject(AliAODMCParticle::StdBranchName());
1982 AliInfo("No MC particle branch found");
1987 AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1988 FindListObject(AliAODMCHeader::StdBranchName());
1989 Float_t vzMC = aodMCheader->GetVtxZ();
1990 if(TMath::Abs(vzMC)>fVertexZCut) return false;
1996 if(recVertex==true){
1997 //if(fAODEvent->IsPileupFromSPD(3,0.8)) return false;
1999 AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
2000 if(!vertex) return false;
2001 TString vtxTitle(vertex->GetTitle());// only allow vertex from tracks, no vertexer z
2002 // Printf("vtxTitle: %s",vtxTitle.Data());
2003 //if (!(vtxTitle.Contains("VertexerTracksWithConstraint"))) return false;
2004 fVertexCheck->Fill(vertex->GetNContributors());
2005 if(vertex->GetNContributors()<=0) return false;
2006 Double_t vzAOD=vertex->GetZ();
2007 // if(TMath::Abs(vzAOD)<1e-9) return false;
2008 if(TMath::Abs(vzAOD)>fVertexZCut) return false;
2010 AliAODVertex* vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();
2011 if(!vertexSPD) return false;
2012 if(vertexSPD->GetNContributors()<=0) return false;
2013 Double_t vzSPD=vertexSPD->GetZ();
2014 //if(TMath::Abs(vzSPD)<1e-9) return false;
2015 if(TMath::Abs(vzSPD)>fVertexZCut) return false;
2017 //check TPC reconstruction: check for corrupted chunks
2018 //better: check TPCvertex, but this is not available yet in AODs
2019 Int_t nAcceptedTracksTPC=0;
2020 Int_t nAcceptedTracksITSTPC=0;
2021 for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2022 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
2023 if (!track) continue;
2024 if(track->TestFilterBit(128) && TMath::Abs(track->Eta())<fEtaCut &&
2025 track->Pt()>fPtMin && track->Pt()<fPtMax)
2026 nAcceptedTracksTPC++;
2027 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut &&
2028 track->Pt()>fPtMin && track->Pt()<fPtMax)
2029 nAcceptedTracksITSTPC++;
2031 fCorruptedChunks->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2033 if(nAcceptedTracksTPC>fNTPC && nAcceptedTracksITSTPC==0)
2034 return false;//most likely corrupted chunk. No ITS tracks are reconstructed
2036 fCorruptedChunksAfter->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2038 //control histograms=================
2040 Int_t ntrackletsAccept=0;
2041 AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
2042 for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
2043 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 &&
2044 TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut) ++ntrackletsAccept;
2046 Int_t nAcceptedTracks=0;
2047 for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2048 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
2049 if (!track) continue;
2050 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
2051 && track->Pt()>fPtMin && track->Pt()<fPtMax) nAcceptedTracks++;
2053 fNContrNtracklets->Fill(ntrackletsAccept,vertexSPD->GetNContributors());
2054 fNContrNtracks->Fill(nAcceptedTracks,vertexSPD->GetNContributors());
2055 //====================================
2062 Printf("Wrong mode!");
2068 //_____________________________________________________________________________
2069 const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins,
2070 const Double_t xmin,
2071 const Double_t xmax)
2073 // returns pointer to an array with can be used to build a logarithmic axis
2074 // it is user responsibility to delete the array
2076 Double_t logxmin = TMath::Log10(xmin);
2077 Double_t logxmax = TMath::Log10(xmax);
2078 Double_t binwidth = (logxmax-logxmin)/nbins;
2080 Double_t *xbins = new Double_t[nbins+1];
2083 for (Int_t i=1;i<=nbins;i++) {
2084 xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth);
2090 //_____________________________________________________________________________
2091 Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis,
2092 const Short_t chargeOthers)
2094 // compute if charge of two particles/tracks has same sign or different sign
2096 if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers);
2098 if(chargeEventAxis<0){
2102 else if(chargeOthers>0){
2107 else if(chargeEventAxis>0){
2111 else if(chargeOthers<0){
2117 AliInfo("Error: Charge not lower nor higher as zero");
2121 AliInfo("Error: Check values of Charge");