11 #include "AliVEvent.h"
12 #include "AliVParticle.h"
14 #include "AliESDEvent.h"
15 #include "AliESDtrack.h"
16 #include "AliMultiplicity.h"
17 #include "AliESDtrackCuts.h"
19 #include "AliAODEvent.h"
20 #include "AliAODTrack.h"
21 #include "AliAODMCParticle.h"
24 #include "AliMCEvent.h"
25 #include "AliMCParticle.h"
26 #include "AliGenEventHeader.h"
28 #include "AliAnalysisManager.h"
30 #include "AliAnalysisTaskMinijet.h"
32 // Analysis task for two-particle correlations using all particles over pt threshold
33 // pt_trig threshold for trigger particle (event axis) and pt_assoc for possible associated particles.
34 // Extract mini-jet yield and fragmentation properties via Delta-Phi histograms of these correlations
35 // post processing of analysis output via macro plot3and2Gaus.C
36 // Can use ESD or AOD, reconstructed and Monte Carlo data as input
37 // Author: Eva Sicking
40 ClassImp(AliAnalysisTaskMinijet)
42 //________________________________________________________________________
43 AliAnalysisTaskMinijet::AliAnalysisTaskMinijet(const char *name)
44 : AliAnalysisTaskSE(name),
57 fSelectParticlesAssoc(1),
71 for(Int_t i = 0;i< 4;i++){
91 fDPhiDEtaEventAxis[i] = 0;
92 fDPhiDEtaEventAxisSeeds[i]= 0;
94 fTriggerNchSeeds[i] = 0;
95 fTriggerTracklet[i] = 0;
99 fNch07Tracklet[i] = 0;
101 fPNch07Tracklet[i] = 0;
102 fDPhiEventAxis[i] = 0;
103 for(Int_t j=0;j<150;j++){
104 fDPhiEventAxisNchBin[i][j] = 0;
105 fDPhiEventAxisNchBinTrig[i][j] = 0;
107 fDPhiEventAxisTrackletBin[i][j] = 0;
108 fDPhiEventAxisTrackletBinTrig[i][j] = 0;
111 DefineOutput(1, TList::Class());
114 //________________________________________________________________________
115 AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet()
119 if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists;
122 //________________________________________________________________________
123 void AliAnalysisTaskMinijet::UserCreateOutputObjects()
127 if(fDebug) Printf("In User Create Output Objects.");
129 fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 150, 0.1, 3.1);
130 fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
131 fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
132 fHistPt->SetMarkerStyle(kFullCircle);
136 fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 150, 0.1, 3.1);
137 fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)");
138 fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
139 fHistPtMC->SetMarkerStyle(kFullCircle);
141 fNmcNch = new TH2F("fNmcNch", "fNmcNch", 100,-0.5,99.5,100,-0.5,99.5);
142 fPNmcNch = new TProfile("pNmcNch", "pNmcNch", 100,-0.5,99.5);
146 fChargedPi0 = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5);
148 TString labels[4]={"ESDrec", "ESDmc", "AODrec", "AODmc"};
150 // for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){
151 for(Int_t i=0;i<4;i++){
153 fVertexZ[i] = new TH1F(Form("fVertexZ%s",labels[i].Data()),
154 Form("fVertexZ%s",labels[i].Data()) ,
156 fPt[i] = new TH1F(Form("fPt%s",labels[i].Data()),
157 Form("fPt%s",labels[i].Data()) ,
159 fEta[i] = new TH1F (Form("fEta%s",labels[i].Data()),
160 Form("fEta%s",labels[i].Data()) ,
162 fPhi[i] = new TH1F(Form("fPhi%s",labels[i].Data()),
163 Form("fPhi%s",labels[i].Data()) ,
164 360, 0.,2*TMath::Pi());
165 fDcaXY[i] = new TH1F(Form("fDcaXY%s",labels[i].Data()),
166 Form("fDcaXY%s",labels[i].Data()) ,
168 fDcaZ[i] = new TH1F(Form("fDcaZ%s",labels[i].Data()),
169 Form("fDcaZ%s",labels[i].Data()) ,
172 fPtSeed[i] = new TH1F(Form("fPSeedt%s",labels[i].Data()),
173 Form("fPtSeed%s",labels[i].Data()) ,
175 fEtaSeed[i] = new TH1F (Form("fEtaSeed%s",labels[i].Data()),
176 Form("fEtaSeed%s",labels[i].Data()) ,
178 fPhiSeed[i] = new TH1F(Form("fPhiSeed%s",labels[i].Data()),
179 Form("fPhiSeed%s",labels[i].Data()) ,
180 360, 0.,2*TMath::Pi());
182 fPtOthers[i] = new TH1F(Form("fPOtherst%s",labels[i].Data()),
183 Form("fPtOthers%s",labels[i].Data()) ,
185 fEtaOthers[i] = new TH1F (Form("fEtaOthers%s",labels[i].Data()),
186 Form("fEtaOthers%s",labels[i].Data()) ,
188 fPhiOthers[i] = new TH1F(Form("fPhiOthers%s",labels[i].Data()),
189 Form("fPhiOthers%s",labels[i].Data()) ,
190 360, 0.,2*TMath::Pi());
191 fPtEtaOthers[i] = new TH2F(Form("fPtEtaOthers%s",labels[i].Data()),
192 Form("fPtEtaOthers%s",labels[i].Data()) ,
193 500, 0., 50, 100, -1., 1);
195 fPhiEta[i] = new TH2F(Form("fPhiEta%s",labels[i].Data()),
196 Form("fPhiEta%s",labels[i].Data()) ,
197 180, 0., 2*TMath::Pi(), 100, -1.,1.);
198 fDPhiDEtaEventAxis[i] = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()),
199 Form("fDPhiDEtaEventAxis%s",labels[i].Data()) ,
200 180, -0.5* TMath::Pi(), 1.5*TMath::Pi(), 200, -4.,4.);
201 fTriggerNch[i] = new TH1F(Form("fTriggerNch%s",labels[i].Data()),
202 Form("fTriggerNch%s",labels[i].Data()) ,
204 fTriggerNchSeeds[i] = new TH2F(Form("fTriggerNchSeeds%s",labels[i].Data()),
205 Form("fTriggerNchSeeds%s",labels[i].Data()) ,
206 250, -0.5, 249.5, 100, -0.5, 99.5);
207 fTriggerTracklet[i] = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()),
208 Form("fTriggerTracklet%s",labels[i].Data()) ,
210 fNch07Nch[i] = new TH2F(Form("fNch07Nch%s",labels[i].Data()),
211 Form("fNch07Nch%s",labels[i].Data()) ,
212 250, -2.5, 247.5,250, -2.5, 247.5);
213 fPNch07Nch[i] = new TProfile(Form("pNch07Nch%s",labels[i].Data()),
214 Form("pNch07Nch%s",labels[i].Data()) ,
216 fNch07Tracklet[i] = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()),
217 Form("fNch07Tracklet%s",labels[i].Data()) ,
218 250, -2.5, 247.5,250, -2.5, 247.5);
219 fNchTracklet[i] = new TH2F(Form("fNchTracklet%s",labels[i].Data()),
220 Form("fNchTracklet%s",labels[i].Data()) ,
221 250, -2.5, 247.5,250, -2.5, 247.5);
222 fPNch07Tracklet[i] = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()),
223 Form("pNch07Tracklet%s",labels[i].Data()) ,
225 fDPhiEventAxis[i] = new TH1F(Form("fDPhiEventAxis%s",
227 Form("fDPhiEventAxis%s",
229 180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
231 for(Int_t j=0;j<150;j++){
232 fDPhiEventAxisNchBin[i][j] = new TH1F(Form("fDPhiEventAxisNchBin%s%02d",
234 Form("fDPhiEventAxisNchBin%s%02d",
235 labels[i].Data(),j) ,
236 180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
237 fDPhiEventAxisNchBinTrig[i][j] = new TH1F(Form("fDPhiEventAxisNchBinTrig%s%02d",
239 Form("fDPhiEventAxisNchBinTrig%s%02d",
240 labels[i].Data(),j) ,
241 180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
243 fDPhiEventAxisTrackletBin[i][j] = new TH1F(Form("fDPhiEventAxisTrackletBin%s%02d",
245 Form("fDPhiEventAxisTrackletBin%s%02d",
246 labels[i].Data(),j) ,
247 180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
249 fDPhiEventAxisTrackletBinTrig[i][j] = new TH1F(Form("fDPhiEventAxisTrackletBinTrig%s%02d",
251 Form("fDPhiEventAxisTrackletBinTrig%s%02d",
252 labels[i].Data(),j) ,
253 180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
257 fHists = new TList();
260 fHists->Add(fHistPt);
262 fHists->Add(fHistPtMC);
263 fHists->Add(fNmcNch);
264 fHists->Add(fPNmcNch);
266 fHists->Add(fChargedPi0);
268 //for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){
269 for(Int_t i=0;i<4;i++){
270 fHists->Add(fVertexZ[i]);
272 fHists->Add(fEta[i]);
273 fHists->Add(fPhi[i]);
274 fHists->Add(fDcaXY[i]);
275 fHists->Add(fDcaZ[i]);
276 fHists->Add(fPtSeed[i]);
277 fHists->Add(fEtaSeed[i]);
278 fHists->Add(fPhiSeed[i]);
279 fHists->Add(fPtOthers[i]);
280 fHists->Add(fEtaOthers[i]);
281 fHists->Add(fPhiOthers[i]);
282 fHists->Add(fPtEtaOthers[i]);
283 fHists->Add(fPhiEta[i]);
284 fHists->Add(fDPhiDEtaEventAxis[i]);
285 fHists->Add(fTriggerNch[i]);
286 fHists->Add(fTriggerNchSeeds[i]);
287 fHists->Add(fTriggerTracklet[i]);
288 fHists->Add(fNch07Nch[i]);
289 fHists->Add(fPNch07Nch[i]);
290 fHists->Add(fNch07Tracklet[i]);
291 fHists->Add(fNchTracklet[i]);
292 fHists->Add(fPNch07Tracklet[i]);
293 fHists->Add(fDPhiEventAxis[i]);
294 for(Int_t j=0;j<150;j++){
295 fHists->Add(fDPhiEventAxisNchBin[i][j]);
296 fHists->Add(fDPhiEventAxisNchBinTrig[i][j]);
298 fHists->Add(fDPhiEventAxisTrackletBin[i][j]);
299 fHists->Add(fDPhiEventAxisTrackletBinTrig[i][j]);
306 //________________________________________________________________________
307 void AliAnalysisTaskMinijet::UserExec(Option_t *)
310 // Called for each event
312 //Printf("Event starts");
314 AliVEvent *event = InputEvent();
316 Error("UserExec", "Could not retrieve event");
320 //get events, either ESD or AOD
321 fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
322 fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
325 //arrays for track properties
330 //number of accepted tracks and tracklets
331 Int_t ntracks = 0; //return value for reading functions for ESD and AOD
332 Int_t *nTracksTracklets = 0; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc,
333 //for real nall=ncharged)
336 //read data and analyse. Implemented for ESD, ESDmc, AOD, AODmc
337 //-------------------------------------------------------------
338 if (fESDEvent && fMode==0) {//ESDs
339 //ESD MC reading and analysis
342 ntracks = LoopESDMC(&pt, &eta, &phi, &charge, &nTracksTracklets); //read mc particles
344 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 1);//analyse
346 if(pt || eta || phi || charge || nTracksTracklets)
347 CleanArrays(pt, eta, phi, charge, nTracksTracklets);// clean up array memory
351 //ESD reading and analysis
354 ntracks = LoopESD(&pt, &eta, &phi,&charge, &nTracksTracklets); //read tracks
357 Printf("ntracks=%d", nTracksTracklets[0]);
358 Printf("ntracklets=%d", nTracksTracklets[1]);
360 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 0); //analyse tracks
362 if(pt || eta || phi || charge || nTracksTracklets)
363 CleanArrays(pt, eta, phi, charge, nTracksTracklets); // clean up array memory
368 if (fAODEvent && fMode==1) {//AODs
370 //AOD MCreading and analysis
373 ntracks = LoopAODMC(&pt, &eta, &phi, &charge, &nTracksTracklets);//read tracks
375 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse
377 if(pt || eta || phi || charge || nTracksTracklets)
378 CleanArrays(pt, eta, phi, charge, nTracksTracklets);// clean up array memory
382 //AOD reading and analysis
385 ntracks = LoopAOD(&pt, &eta, &phi, &charge, &nTracksTracklets);//read tracks
388 Printf("ntracks=%d", nTracksTracklets[0]);
389 Printf("ntracklets=%d", nTracksTracklets[1]);
391 Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 2);//analyse
393 if(pt || eta || phi || charge || nTracksTracklets)
394 CleanArrays(pt, eta, phi, charge, nTracksTracklets);// clean up array memory
398 //-------------------------------------------------------------
402 // PostData(1, fHists);
407 //________________________________________________________________________
408 Int_t AliAnalysisTaskMinijet::LoopESD(Float_t **ptArray, Float_t ** etaArray,
409 Float_t ** phiArray, Short_t ** chargeArray,
410 Int_t **nTracksTracklets )
412 // gives back the number of esd tracks and pointer to arrays with track
413 // properties (pt, eta, phi)
414 // Uses TPC tracks with SPD vertex from now on
416 // Retreive the number of all tracks for this event.
417 Int_t ntracks = fESDEvent->GetNumberOfTracks();
418 if(fDebug)Printf("all ESD tracks: %d", ntracks);
420 const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
422 //first loop to check how many tracks are accepted
424 Int_t nAcceptedTracks=0;
425 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
427 AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
429 Error("UserExec", "Could not receive track %d", iTracks);
432 //if(!fCuts->AcceptTrack(esdTrack)) continue;
435 // use TPC only tracks with non default SPD vertex
436 AliESDtrack *track = AliESDtrackCuts::
437 GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
439 if(!fCuts->AcceptTrack(track)) continue;// apply TPC track cuts before constrain to SPD vertex
441 // only constrain tracks above threshold
442 AliExternalTrackParam exParam;
443 // take the B-field from the ESD, no 3D fieldMap available at this point
444 Bool_t relate = false;
445 relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
451 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
452 exParam.GetCovariance());
454 else continue;// only if tracks have pt<=0
457 if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.2 && track->Pt()<200.)
460 // TPC only track memory needs to be cleaned up
468 *ptArray = new Float_t[nAcceptedTracks];
469 *etaArray = new Float_t[nAcceptedTracks];
470 *phiArray = new Float_t[nAcceptedTracks];
471 *chargeArray = new Short_t[nAcceptedTracks];
472 *nTracksTracklets = new Int_t[3]; //ntracksAccepted, ntracklets
474 //check if event is pile up or no tracks are accepted, return to user exec
475 if(fESDEvent->IsPileupFromSPD(3,0,8)) return 0;
476 if(nAcceptedTracks==0) return 0;
478 //accept event only, if vertex is good and is within fVertexZcut region
479 const AliESDVertex* vertexESD = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
481 if(!vertexESD) return 0;
482 if(vertexESD->GetNContributors()<=0)return 0;
483 Float_t fVz= vertexESD->GetZ();
484 if(TMath::Abs(fVz)>fVertexZCut) return 0;
485 fVertexZ[0]->Fill(fVz);
487 //variables for DCA QA check per track
492 Int_t iAcceptedTrack=0;
493 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
494 AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
496 Error("UserExec", "Could not receive track %d", iTracks);
499 //if(!fCuts->AcceptTrack(esdTrack)) continue;
501 // create a tpc only track with SPD vertex
502 AliESDtrack *track = AliESDtrackCuts::
503 GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
505 if(!fCuts->AcceptTrack(track)) continue; // apply TPC track cuts before constrain to SPD vertex
508 // only constrain tracks above threshold
509 AliExternalTrackParam exParam;
510 // take the B-field from the ESD, no 3D fieldMap available at this point
511 Bool_t relate = false;
512 relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
518 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
519 exParam.GetCovariance());
524 if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.2 && track->Pt()<200.){
525 (*ptArray)[iAcceptedTrack] = track->Pt();
526 (*etaArray)[iAcceptedTrack] = track->Eta();
527 (*phiArray)[iAcceptedTrack] = track->Phi();
528 (*chargeArray)[iAcceptedTrack++] = track->Charge();
529 fHistPt->Fill(track->Pt());
533 track->GetImpactParameters(fXY,fZ);
534 fDcaXY[0]->Fill(fXY);
539 // TPC only track memory needs to be cleaned up
548 Int_t ntrackletsAccept=0;
549 AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
550 Int_t ntracklets = mult->GetNumberOfTracklets();
551 for(Int_t i=0;i< ntracklets;i++){
552 if(mult->GetDeltaPhi(i)<0.05){
557 (*nTracksTracklets)[0] = nAcceptedTracks;
558 (*nTracksTracklets)[1] = ntrackletsAccept;
559 (*nTracksTracklets)[2] = nAcceptedTracks;//in order to be compatible with mc analysis
560 //where here also neutral particles are counted.
564 //Printf("Number of MC particles from ESDMC = %d",fNMcPrimAccept);
565 //Printf("Number of tracks from ESD = %d",nAcceptedTracks);
566 fNmcNch->Fill(fNMcPrimAccept,nAcceptedTracks);
567 fPNmcNch->Fill(fNMcPrimAccept,nAcceptedTracks);
568 return fNMcPrimAccept; // also possible to use reconstructed Nch -> return nAcceptedTracks;
571 fVzEvent=fVz; // needed for correction maps
572 return nAcceptedTracks;
576 //________________________________________________________________________
577 Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray,
578 Float_t ** phiArray, Short_t ** chargeArray,
579 Int_t ** nTracksTracklets)
581 // gives back the number of charged prim MC particle and pointer to arrays
582 // with particle properties (pt, eta, phi)
584 // Printf("Event starts: Loop ESD MC");
586 AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
588 Error("UserExec", "Could not retrieve MC event");
592 AliStack* stack = MCEvent()->Stack();
595 Int_t ntracks = mcEvent->GetNumberOfTracks();
596 if(fDebug)Printf("MC particles: %d", ntracks);
599 //----------------------------------
600 //first track loop to check how many chared primary tracks are available
601 Int_t nChargedPrimaries=0;
602 Int_t nAllPrimaries=0;
604 Int_t nPseudoTracklets=0;
605 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
606 AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
608 Error("UserExec", "Could not receive track %d", iTracks);
613 if(//count also charged particles in case of fSelectParticles==2 (only neutral)
614 !SelectParticlePlusCharged(
617 stack->IsPhysicalPrimary(track->Label())
623 // Printf("fSelectParticles= %d\n", fSelectParticles);
624 //count number of pseudo tracklets
625 if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.035) ++nPseudoTracklets;
626 //same cuts as on ESDtracks
627 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue;
629 //count all primaries
631 //count charged primaries
632 if (track->Charge() != 0) ++nChargedPrimaries;
634 // Printf("PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
638 //----------------------------------
640 // Printf("All in acceptance=%d",nAllPrimaries);
641 // Printf("Charged in acceptance =%d",nChargedPrimaries);
643 fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
645 if(fSelectParticles!=2){
646 *ptArray = new Float_t[nAllPrimaries];
647 *etaArray = new Float_t[nAllPrimaries];
648 *phiArray = new Float_t[nAllPrimaries];
649 *chargeArray = new Short_t[nAllPrimaries];
652 *ptArray = new Float_t[nAllPrimaries-nChargedPrimaries];
653 *etaArray = new Float_t[nAllPrimaries-nChargedPrimaries];
654 *phiArray = new Float_t[nAllPrimaries-nChargedPrimaries];
655 *chargeArray = new Short_t[nAllPrimaries-nChargedPrimaries];
658 *nTracksTracklets = new Int_t[3];
660 if(nAllPrimaries==0) return 0;
661 if(nChargedPrimaries==0) return 0;
665 AliGenEventHeader* header = MCEvent()->GenEventHeader();
667 header->PrimaryVertex(mcV);
668 if(TMath::Abs(mcV[0])<1e-8 && TMath::Abs(mcV[0])<1e-8 && TMath::Abs(mcV[0])<1e-8) return 0;
669 Float_t vzMC = mcV[2];
670 if(TMath::Abs(vzMC)>fVertexZCut) return 0;
671 fVertexZ[1]->Fill(vzMC);
676 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
677 AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
679 Error("UserExec", "Could not receive track %d", iTracks);
686 stack->IsPhysicalPrimary(track->Label())
692 //same cuts as on ESDtracks
693 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2
694 || track->Pt()>200) continue;
696 // Printf("After: PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
699 fHistPtMC->Fill(track->Pt());
700 //fills arrays with track properties
701 (*ptArray)[iChargedPiK] = track->Pt();
702 (*etaArray)[iChargedPiK] = track->Eta();
703 (*phiArray)[iChargedPiK] = track->Phi();
704 (*chargeArray)[iChargedPiK++] = track->Charge();
708 (*nTracksTracklets)[0] = nChargedPrimaries;
709 (*nTracksTracklets)[1] = nPseudoTracklets;
710 if(fSelectParticles!=2){
711 (*nTracksTracklets)[2] = nAllPrimaries;
714 (*nTracksTracklets)[2] = nAllPrimaries-nChargedPrimaries; // only neutral
717 // Printf("Number of MC particles = %d",nChargedPrimaries);
718 fNMcPrimAccept=nChargedPrimaries;
719 return nChargedPrimaries;
723 //________________________________________________________________________
724 Int_t AliAnalysisTaskMinijet::LoopAOD(Float_t **ptArray, Float_t ** etaArray,
725 Float_t ** phiArray, Short_t ** chargeArray,
726 Int_t ** nTracksTracklets)
728 // gives back the number of AOD tracks and pointer to arrays with track
729 // properties (pt, eta, phi)
732 // Retreive the number of tracks for this event.
733 Int_t ntracks = fAODEvent->GetNumberOfTracks();
734 if(fDebug) Printf("AOD tracks: %d", ntracks);
737 Int_t nAcceptedTracks=0;
738 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
739 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
741 Error("UserExec", "Could not receive track %d", iTracks);
744 // filter bit needs to be ajusted to TPC only tracks with SPD vertex as soon as AOD are available (so far ESD track cuts ITS TPC 2010)
745 if(track->TestFilterBit(16) && TMath::Abs(track->Eta())<fEtaCut
746 && track->Pt()>0.2 && track->Pt()<200.) nAcceptedTracks++;
749 *ptArray = new Float_t[nAcceptedTracks];
750 *etaArray = new Float_t[nAcceptedTracks];
751 *phiArray = new Float_t[nAcceptedTracks];
752 *chargeArray = new Short_t[nAcceptedTracks];
753 *nTracksTracklets = new Int_t[3]; //here, third entry only copy of first (compatibility for MC)
756 if(nAcceptedTracks==0) return 0;
757 AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
759 // TODO: need to check how to implement IsPileupFromSPD(3,0.8)
760 // function of esd event
761 // first solution: exclude this call in esd loop for comparison (QA)
763 if(!vertex) return 0;
764 Double_t vzAOD=vertex->GetZ();
765 if(vertex->GetNContributors()<=0) return 0;
766 if(TMath::Abs(vzAOD)<1e-9) return 0;
768 if(TMath::Abs(vzAOD)>fVertexZCut) return 0;
769 fVertexZ[2]->Fill(vzAOD);
771 // Track loop to fill a pT spectrum
772 Int_t iAcceptedTracks=0;
773 for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
774 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
776 Error("UserExec", "Could not receive track %d", iTracks);
779 if(!track->TestFilterBit(16) || TMath::Abs(track->Eta())>fEtaCut
780 || track->Pt()<0.2 || track->Pt()>200.) continue;
781 fHistPt->Fill(track->Pt());
783 //fills arrays with track properties
784 (*ptArray)[iAcceptedTracks] = track->Pt();
785 (*etaArray)[iAcceptedTracks] = track->Eta();
786 (*phiArray)[iAcceptedTracks] = track->Phi();
787 (*chargeArray)[iAcceptedTracks++] = track->Charge();
792 Int_t ntrackletsAccept=0;
793 AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
794 for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
795 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05){
800 (*nTracksTracklets)[0] = nAcceptedTracks;
801 (*nTracksTracklets)[1] = ntrackletsAccept;
802 (*nTracksTracklets)[2] = nAcceptedTracks;
805 fNmcNch->Fill(fNMcPrimAccept,nAcceptedTracks);
806 return fNMcPrimAccept;
809 return nAcceptedTracks;
813 //________________________________________________________________________
814 Int_t AliAnalysisTaskMinijet::LoopAODMC(Float_t **ptArray, Float_t ** etaArray,
815 Float_t ** phiArray, Short_t ** chargeArray,
816 Int_t ** nTracksTracklets)
818 // gives back the number of AOD MC particles and pointer to arrays with particle
819 // properties (pt, eta, phi)
821 //retreive MC particles from event
822 TClonesArray *mcArray = (TClonesArray*)fAODEvent->
823 FindListObject(AliAODMCParticle::StdBranchName());
825 Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
829 Int_t ntracks = mcArray->GetEntriesFast();
830 if(fDebug)Printf("MC particles: %d", ntracks);
833 // Track loop: chek how many particles will be accepted
835 Int_t nChargedPrim=0;
837 Int_t nPseudoTracklets=0;
838 for (Int_t it = 0; it < ntracks; it++) {
839 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
841 Error("UserExec", "Could not receive track %d", it);
845 if(!SelectParticlePlusCharged(
848 track->IsPhysicalPrimary()
853 if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.035)++nPseudoTracklets;
854 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue;
855 // if(TMath::Abs(track->Eta())<1e-9 && TMath::Abs(track->Phi())<1e-9)continue;
857 //same cuts as in ESD filter
858 //if(!track->TestBit(16))continue; //cuts set in ESD filter during AOD generation
861 if(track->Charge()!=0) nChargedPrim++;
863 // Printf("eta=%f,phi=%f,pt=%f",track->Eta(),track->Phi(),track->Pt());
864 // Printf("prim=%d,pdg=%d,charge=%d",track->IsPhysicalPrimary(),track->PdgCode(),track->Charge());
867 if(nAllPrim==1) vzMC= track->Zv(); // check only one time. (only one vertex per event allowed)
870 //generate array with size of number of accepted tracks
871 fChargedPi0->Fill(nAllPrim,nChargedPrim);
873 if(fSelectParticles!=2){
874 *ptArray = new Float_t[nAllPrim];
875 *etaArray = new Float_t[nAllPrim];
876 *phiArray = new Float_t[nAllPrim];
877 *chargeArray = new Short_t[nAllPrim];
880 *ptArray = new Float_t[nAllPrim-nChargedPrim];
881 *etaArray = new Float_t[nAllPrim-nChargedPrim];
882 *phiArray = new Float_t[nAllPrim-nChargedPrim];
883 *chargeArray = new Short_t[nAllPrim-nChargedPrim];
887 *nTracksTracklets = new Int_t[3];
889 //Printf("nAllPrim=%d", nAllPrim);
890 //Printf("nChargedPrim=%d", nChargedPrim);
892 if(nAllPrim==0) return 0;
893 if(nChargedPrim==0) return 0;
896 if(TMath::Abs(vzMC)>fVertexZCut) return 0;
897 fVertexZ[3]->Fill(vzMC);
900 // Track loop: fill arrays for accepted tracks
901 Int_t iChargedPrim=0;
902 for (Int_t it = 0; it < ntracks; it++) {
903 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
905 Error("UserExec", "Could not receive track %d", it);
912 track->IsPhysicalPrimary()
918 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue;
919 //if(TMath::Abs(track->Eta())<1e-8 && TMath::Abs(track->Phi())<1e-8)continue;
921 //Printf("eta=%f,phi=%f,pt=%f",track->Eta(),track->Phi(),track->Pt());
922 //Printf("prim=%d,pdg=%d,charge=%d",track->IsPhysicalPrimary(),track->PdgCode(),track->Charge());
924 //if(track->TestBit(16))continue;
926 fHistPtMC->Fill(track->Pt());
927 (*ptArray)[iChargedPrim] = track->Pt();
928 (*etaArray)[iChargedPrim] = track->Eta();
929 (*phiArray)[iChargedPrim] = track->Phi();
930 (*chargeArray)[iChargedPrim++] = track->Charge();
934 (*nTracksTracklets)[0] = nChargedPrim;
935 (*nTracksTracklets)[1] = nPseudoTracklets;
936 if(fSelectParticles!=2){
937 (*nTracksTracklets)[2] = nAllPrim;
940 (*nTracksTracklets)[2] = nAllPrim-nChargedPrim; // only neutral
943 fNMcPrimAccept=nChargedPrim;
945 // Printf("ntracks=%d", nChargedPrim);
949 //________________________________________________________________________
950 void AliAnalysisTaskMinijet::Analyse(const Float_t *pt, const Float_t *eta, const Float_t *phi,
951 const Short_t *charge, Int_t ntracksCharged,
952 Int_t ntracklets, const Int_t nAll, Int_t mode)
955 // analyse track properties (comming from either ESDs or AODs) in order to compute
956 // mini jet activity (chared tracks) as function of charged multiplicity
958 // ntracks and ntracklets are already the number of accepted tracks and tracklets
961 Printf("In Analysis\n");
962 Printf("nAll=%d",nAll);
963 Printf("nCharged=%d",ntracksCharged);
966 Float_t ptEventAxis=0; // pt leading
967 Float_t etaEventAxis=0; // eta leading
968 Float_t phiEventAxis=0; // phi leading
970 Float_t ptOthers = 0; // pt others // for all other tracks around event axis -> see loop
971 Float_t etaOthers = 0; // eta others
972 Float_t phiOthers = 0; // phi others
974 Int_t *pindexInnerEta = new Int_t[nAll+1]; // Fix for coverity defect complaining in TMath:Sort in line 1022
975 Float_t *ptInnerEta = new Float_t[nAll+1]; // ToDo: Proper fix needs to be investigated
977 for (Int_t i = 0; i < nAll; i++) {
978 //filling of simple check plots
979 fPt[mode] -> Fill( pt[i]);
980 fEta[mode] -> Fill(eta[i]);
981 fPhi[mode] -> Fill(phi[i]);
982 fPhiEta[mode]-> Fill(phi[i], eta[i]);
984 pindexInnerEta[i]=0; //set all values to zero
985 //fill new array for eta check
986 ptInnerEta[i]= pt[i];
991 // define event axis: leading or random track (with pt>fTriggerPtCut)
992 // ---------------------------------------
993 Int_t highPtTracks=0;
994 Int_t highPtTracksInnerEta=0;
997 //count high pt tracks and high pt tracks in acceptance for seeds
998 for(Int_t i=0;i<nAll;i++){
1000 if(TMath::Abs(eta[i])<0.9){
1004 if(pt[i]>fTriggerPtCut) {
1008 // seed should be place in middle of acceptance, that complete cone is in acceptance
1009 if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
1011 // Printf("eta=%f", eta[i]);
1012 highPtTracksInnerEta++;
1020 //sort array in order to get highest pt tracks first
1021 //index can be computed with pindexInnerEta[number]
1022 if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1025 // plot of multiplicity distributions
1026 fNch07Nch[mode]->Fill(ntracksCharged, highPtTracksInnerEta);
1027 fPNch07Nch[mode]->Fill(ntracksCharged, highPtTracksInnerEta);
1029 fNch07Tracklet[mode]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1030 fNchTracklet[mode]->Fill(ntracklets, ntracksCharged);
1031 fPNch07Tracklet[mode]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1034 //analysis can only be performed with event axis, defined by high pt track
1037 if(highPtTracks>0 && highPtTracksInnerEta>0){
1039 //Printf("%s:%d",(char*)__FILE__,__LINE__);
1040 //check setter of event axis
1041 //default option: random=1,
1042 //second option:leading=0
1045 // if(fLeadingOrRandom==0)axis=0;
1046 // else if (fLeadingOrRandom==1)axis= Int_t((highPtTracksInnerEta)*gRandom->Rndm());
1047 // else Printf("Wrong settings for event axis.");
1050 // Printf("Axis tracks has pT=%f, test=%f", ptInnerEta[pindexInnerEta[axis]], pt[pindexInnerEta[axis]]);
1051 // Printf("Axis tracks has eta=%f", eta[pindexInnerEta[axis]]);
1054 //---------------------------------------
1055 //Printf("Number of seeds = %d",highPtTracksInnerEta);
1058 // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
1059 for(Int_t axis=0;axis<highPtTracksInnerEta; axis++){
1061 if(ntracksCharged>1){ // require at least two tracks (leading and prob. accosicates)
1063 //EventAxisRandom track properties
1064 ptEventAxis = pt [pindexInnerEta[axis]];
1065 etaEventAxis = eta[pindexInnerEta[axis]];
1066 phiEventAxis = phi[pindexInnerEta[axis]];
1067 fPtSeed[mode] -> Fill( ptEventAxis);
1068 fEtaSeed[mode] -> Fill(etaEventAxis);
1069 fPhiSeed[mode] -> Fill(phiEventAxis);
1071 //track loop for event propoerties around event axis with pt>triggerPtCut
1072 //loop only over already accepted tracks except event axis
1073 // if(ptEventAxis>fTriggerPtCut && ptEventAxis < 0.8){
1074 if(ptEventAxis>fTriggerPtCut){
1076 for (Int_t iTrack = 0; iTrack < nAll; iTrack++) {
1078 if(pindexInnerEta[axis]==iTrack)continue; // no double counting
1080 if(fSelectParticlesAssoc==1){
1081 if(charge[iTrack]==0)continue;
1083 if(fSelectParticlesAssoc==2){
1084 if(charge[iTrack]!=0)continue;
1086 // Printf("Charge=%d", charge[iTrack]);
1089 ptOthers = pt [iTrack];
1090 etaOthers = eta[iTrack];
1091 phiOthers = phi[iTrack];
1094 //if(ptOthers>0.4 && ptOthers<0.5){ // only tracks which fullfill associate pt cut
1095 if(ptOthers>fAssociatePtCut){ // only tracks which fullfill associate pt cut
1097 //plot only properties of tracks with pt>ptassoc
1098 fPtOthers[mode] -> Fill( ptOthers);
1099 fEtaOthers[mode] -> Fill(etaOthers);
1100 fPhiOthers[mode] -> Fill(phiOthers);
1101 fPtEtaOthers[mode] -> Fill(ptOthers, etaOthers);
1103 Float_t dPhi = phiOthers-phiEventAxis;
1104 if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1105 else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1106 Float_t dEta=etaOthers-etaEventAxis;
1108 fDPhiDEtaEventAxis[mode]->Fill(dPhi, dEta);
1110 fDPhiEventAxis[mode]->Fill(dPhi);
1111 if(ntracksCharged<150){ // in case of MC data, ntracksCharged represents Nmcprim for ESD and MC loop
1112 fDPhiEventAxisNchBin[mode][ntracksCharged]->Fill(dPhi);
1113 if(ptOthers>fTriggerPtCut && TMath::Abs(etaOthers)<fEtaCutSeed)
1114 fDPhiEventAxisNchBinTrig[mode][ntracksCharged]->Fill(dPhi);
1118 fDPhiEventAxisTrackletBin[mode][ntracklets]->Fill(dPhi);
1119 if(ptOthers>fTriggerPtCut && TMath::Abs(etaOthers)<fEtaCutSeed )
1120 fDPhiEventAxisTrackletBinTrig[mode][ntracklets]->Fill(dPhi);
1123 }//tracks fulfill assoc track cut
1125 }// end of inner track loop
1128 // fill histogram with number of tracks (pt>fAssociatePtCut) around event axis
1129 // how often is there a trigger particle at a certain Nch bin
1132 }//if track pt is at least trigger pt
1134 }//loop over all highPtInnerEta tracks
1137 } //if there are more than 1 track
1139 fTriggerNch[mode]->Fill(ntracksCharged,highPtTracksInnerEta);
1140 fTriggerNchSeeds[mode]->Fill(ntracksCharged,highPtTracksInnerEta);
1141 fTriggerTracklet[mode]->Fill(ntracklets);
1143 }//if there is at least one high pt track
1146 if(pindexInnerEta){// clean up array memory used for TMath::Sort
1147 delete[] pindexInnerEta;
1151 if(ptInnerEta){// clean up array memory used for TMath::Sort
1152 delete[] ptInnerEta;
1156 PostData(1, fHists);
1160 //________________________________________________________________________
1161 void AliAnalysisTaskMinijet::Terminate(Option_t*)
1163 //terminate function is called at the end
1164 //can be used to draw histograms etc.
1168 //________________________________________________________________________
1169 Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, const Bool_t prim)
1171 //selection of mc particle
1172 //fSelectParticles=0: use charged primaries and pi0 and k0
1173 //fSelectParticles=1: use only charged primaries
1174 //fSelectParticles=2: use only pi0 and k0
1176 if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1178 if(!(pdg==111||pdg==130||pdg==310))
1187 else if(fSelectParticles==1){
1188 if (charge==0 || !prim){
1194 Printf("Error: wrong selection of charged/pi0/k0");
1202 //________________________________________________________________________
1203 Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
1205 //selection of mc particle
1206 //fSelectParticles=0: use charged primaries and pi0 and k0
1207 //fSelectParticles=1: use only charged primaries
1208 //fSelectParticles=2: use only pi0 and k0
1210 if(fSelectParticles==0){
1212 if(!(pdg==111||pdg==130||pdg==310))
1221 else if (fSelectParticles==1){
1222 if (charge==0 || !prim){
1226 else if(fSelectParticles==2){
1227 if(!(pdg==111||pdg==130||pdg==310))
1235 //________________________________________________________________________
1236 void AliAnalysisTaskMinijet::CleanArrays(const Float_t* pt,
1239 const Short_t * charge,
1240 const Int_t* nTracksTracklets)
1242 //clean up of memory used for arrays of track properties
1260 if(nTracksTracklets){
1261 delete[] nTracksTracklets;