1 /*************************************************************************
4 * Task for Jet Chemistry Analysis in PWG4 Jet Task Force Train *
7 *************************************************************************/
9 /**************************************************************************
10 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
12 * Author: The ALICE Off-line Project. *
13 * Contributors are mentioned in the code where appropriate. *
15 * Permission to use, copy, modify and distribute this software and its *
16 * documentation strictly for non-commercial purposes is hereby granted *
17 * without fee, provided that the above copyright notice appears in all *
18 * copies and that both the copyright notice and this permission notice *
19 * appear in the supporting documentation. The authors make no claims *
20 * about the suitability of this software for any purpose. It is *
21 * provided "as is" without express or implied warranty. *
22 **************************************************************************/
29 #include "THnSparse.h"
31 #include "AliAnalysisHelperJetTasks.h"
32 #include "AliAnalysisManager.h"
33 #include "AliAODHandler.h"
34 #include "AliAODInputHandler.h"
35 #include "AliESDEvent.h"
36 #include "AliGenPythiaEventHeader.h"
37 #include "AliGenHijingEventHeader.h"
39 #include "AliAODEvent.h"
40 #include "AliAODJet.h"
42 #include "AliAODTrack.h"
46 #include "AliExternalTrackParam.h"
48 #include "AliAnalysisTaskJetChem.h"
50 ClassImp(AliAnalysisTaskJetChem)
52 //____________________________________________________________________________
53 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem()
54 : AliAnalysisTaskFragmentationFunction()
59 ,fFFHistosRecCutsK0Evt(0)
60 ,fFFHistosIMK0AllEvt(0)
63 ,fFFHistosPhiCorrIMK0(0)
82 ,fPhiCorrIMNBinsPhi(0)
85 ,fPhiCorrIMNBinsInvM(0)
91 // default constructor
94 //__________________________________________________________________________________________
95 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const char *name)
96 : AliAnalysisTaskFragmentationFunction(name)
101 ,fFFHistosRecCutsK0Evt(0)
102 ,fFFHistosIMK0AllEvt(0)
104 ,fFFHistosIMK0Cone(0)
105 ,fFFHistosPhiCorrIMK0(0)
121 ,fPhiCorrIMNBinsPt(0)
124 ,fPhiCorrIMNBinsPhi(0)
127 ,fPhiCorrIMNBinsInvM(0)
128 ,fPhiCorrIMInvMMin(0)
129 ,fPhiCorrIMInvMMax(0)
135 DefineOutput(1,TList::Class());
138 //__________________________________________________________________________________________________________________________
139 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const AliAnalysisTaskJetChem ©)
140 : AliAnalysisTaskFragmentationFunction()
141 ,fK0Type(copy.fK0Type)
142 ,fFilterMaskK0(copy.fFilterMaskK0)
143 ,fListK0s(copy.fListK0s)
144 ,fV0QAK0(copy.fV0QAK0)
145 ,fFFHistosRecCutsK0Evt(copy.fFFHistosRecCutsK0Evt)
146 ,fFFHistosIMK0AllEvt(copy.fFFHistosIMK0AllEvt)
147 ,fFFHistosIMK0Jet(copy.fFFHistosIMK0Jet)
148 ,fFFHistosIMK0Cone(copy.fFFHistosIMK0Cone)
149 ,fFFHistosPhiCorrIMK0(copy.fFFHistosPhiCorrIMK0)
150 ,fFFIMNBinsJetPt(copy.fFFIMNBinsJetPt)
151 ,fFFIMJetPtMin(copy.fFFIMJetPtMin)
152 ,fFFIMJetPtMax(copy.fFFIMJetPtMax)
153 ,fFFIMNBinsInvM(copy.fFFIMNBinsInvM)
154 ,fFFIMInvMMin(copy.fFFIMInvMMin)
155 ,fFFIMInvMMax(copy.fFFIMInvMMax)
156 ,fFFIMNBinsPt(copy.fFFIMNBinsPt)
157 ,fFFIMPtMin(copy.fFFIMPtMin)
158 ,fFFIMPtMax(copy.fFFIMPtMax)
159 ,fFFIMNBinsXi(copy.fFFIMNBinsXi)
160 ,fFFIMXiMin(copy.fFFIMXiMin)
161 ,fFFIMXiMax(copy.fFFIMXiMax)
162 ,fFFIMNBinsZ(copy.fFFIMNBinsZ)
163 ,fFFIMZMin(copy.fFFIMZMin)
164 ,fFFIMZMax(copy.fFFIMZMax)
165 ,fPhiCorrIMNBinsPt(copy.fPhiCorrIMNBinsPt)
166 ,fPhiCorrIMPtMin(copy.fPhiCorrIMPtMin)
167 ,fPhiCorrIMPtMax(copy.fPhiCorrIMPtMax)
168 ,fPhiCorrIMNBinsPhi(copy.fPhiCorrIMNBinsPhi)
169 ,fPhiCorrIMPhiMin(copy.fPhiCorrIMPhiMin)
170 ,fPhiCorrIMPhiMax(copy.fPhiCorrIMPhiMax)
171 ,fPhiCorrIMNBinsInvM(copy.fPhiCorrIMNBinsInvM)
172 ,fPhiCorrIMInvMMin(copy.fPhiCorrIMInvMMin)
173 ,fPhiCorrIMInvMMax(copy.fPhiCorrIMInvMMax)
174 ,fh1K0Mult(copy.fh1K0Mult)
175 ,fh1dPhiJetK0(copy.fh1dPhiJetK0)
181 // _________________________________________________________________________________________________________________________________
182 AliAnalysisTaskJetChem& AliAnalysisTaskJetChem::operator=(const AliAnalysisTaskJetChem& o)
187 AliAnalysisTaskFragmentationFunction::operator=(o);
190 fFilterMaskK0 = o.fFilterMaskK0;
191 fListK0s = o.fListK0s;
193 fFFHistosRecCutsK0Evt = o.fFFHistosRecCutsK0Evt;
194 fFFHistosIMK0AllEvt = o.fFFHistosIMK0AllEvt;
195 fFFHistosIMK0Jet = o.fFFHistosIMK0Jet;
196 fFFHistosIMK0Cone = o.fFFHistosIMK0Cone;
197 fFFHistosPhiCorrIMK0 = o.fFFHistosPhiCorrIMK0;
198 fFFIMNBinsJetPt = o.fFFIMNBinsJetPt;
199 fFFIMJetPtMin = o.fFFIMJetPtMin;
200 fFFIMJetPtMax = o.fFFIMJetPtMax;
201 fFFIMNBinsPt = o.fFFIMNBinsPt;
202 fFFIMPtMin = o.fFFIMPtMin;
203 fFFIMPtMax = o.fFFIMPtMax;
204 fFFIMNBinsXi = o.fFFIMNBinsXi;
205 fFFIMXiMin = o.fFFIMXiMin;
206 fFFIMXiMax = o.fFFIMXiMax;
207 fFFIMNBinsZ = o.fFFIMNBinsZ;
208 fFFIMZMin = o.fFFIMZMin;
209 fFFIMZMax = o.fFFIMZMax;
210 fPhiCorrIMNBinsPt = o.fPhiCorrIMNBinsPt;
211 fPhiCorrIMPtMin = o.fPhiCorrIMPtMin;
212 fPhiCorrIMPtMax = o.fPhiCorrIMPtMax;
213 fPhiCorrIMNBinsPhi = o.fPhiCorrIMNBinsPhi;
214 fPhiCorrIMPhiMin = o.fPhiCorrIMPhiMin;
215 fPhiCorrIMPhiMax = o.fPhiCorrIMPhiMax;
216 fPhiCorrIMNBinsInvM = o.fPhiCorrIMNBinsInvM;
217 fPhiCorrIMInvMMin = o.fPhiCorrIMInvMMin;
218 fPhiCorrIMInvMMax = o.fPhiCorrIMInvMMax;
219 fh1K0Mult = o.fh1K0Mult;
220 fh1dPhiJetK0 = o.fh1dPhiJetK0;
226 //_______________________________________________
227 AliAnalysisTaskJetChem::~AliAnalysisTaskJetChem()
232 if(fTracksRecCuts) delete fTracksRecCuts;
233 if(fJetsRecCuts) delete fJetsRecCuts;
234 if(fListK0s) delete fListK0s;
235 if(fCommonHistList) delete fCommonHistList;
239 //________________________________________________________________________________________________________________________________
240 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const char* name,
241 Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,
242 Int_t nInvMass, Float_t invMassMin, Float_t invMassMax,
243 Int_t nPt, Float_t ptMin, Float_t ptMax,
244 Int_t nXi, Float_t xiMin, Float_t xiMax,
245 Int_t nZ , Float_t zMin , Float_t zMax )
250 ,fNBinsInvMass(nInvMass)
251 ,fInvMassMin(invMassMin)
252 ,fInvMassMax(invMassMax)
268 // default constructor
272 //______________________________________________________________________________________________________________
273 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const AliFragFuncHistosInvMass& copy)
275 ,fNBinsJetPt(copy.fNBinsJetPt)
276 ,fJetPtMin(copy.fJetPtMin)
277 ,fJetPtMax(copy.fJetPtMax)
278 ,fNBinsInvMass(copy.fNBinsInvMass)
279 ,fInvMassMin(copy.fInvMassMin)
280 ,fInvMassMax(copy.fInvMassMax)
281 ,fNBinsPt(copy.fNBinsPt)
284 ,fNBinsXi(copy.fNBinsXi)
287 ,fNBinsZ(copy.fNBinsZ)
290 ,fh3TrackPt(copy.fh3TrackPt)
293 ,fh1JetPt(copy.fh1JetPt)
294 ,fNameFF(copy.fNameFF)
299 //______________________________________________________________________________________________________________________________________________________________________
300 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& o)
305 TObject::operator=(o);
306 fNBinsJetPt = o.fNBinsJetPt;
307 fJetPtMin = o.fJetPtMin;
308 fJetPtMax = o.fJetPtMax;
309 fNBinsInvMass = o.fNBinsInvMass;
310 fInvMassMin = o.fInvMassMin;
311 fInvMassMax = o.fInvMassMax;
312 fNBinsPt = o.fNBinsPt;
315 fNBinsXi = o.fNBinsXi;
321 fh3TrackPt = o.fh3TrackPt;
324 fh1JetPt = o.fh1JetPt;
331 //___________________________________________________________________________
332 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::~AliFragFuncHistosInvMass()
336 if(fh1JetPt) delete fh1JetPt;
337 if(fh3TrackPt) delete fh3TrackPt;
338 if(fh3Xi) delete fh3Xi;
339 if(fh3Z) delete fh3Z;
342 //_________________________________________________________________
343 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::DefineHistos()
347 fh1JetPt = new TH1F(Form("fh1FFJetPtIM%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
348 fh3TrackPt = new TH3F(Form("fh3FFTrackPtIM%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsPt, fPtMin, fPtMax);
349 fh3Xi = new TH3F(Form("fh3FFXiIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsXi, fXiMin, fXiMax);
350 fh3Z = new TH3F(Form("fh3FFZIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsZ, fZMin, fZMax);
352 AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{t} (GeV/c)", "entries");
353 AliAnalysisTaskJetChem::SetProperties(fh3TrackPt,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","p_{t} (GeV/c)");
354 AliAnalysisTaskJetChem::SetProperties(fh3Xi,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","#xi");
355 AliAnalysisTaskJetChem::SetProperties(fh3Z,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","z");
358 //________________________________________________________________________________________________________________________________
359 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::FillFF(Float_t trackPt, Float_t invM, Float_t jetPt, Bool_t incrementJetPt)
363 if(incrementJetPt) fh1JetPt->Fill(jetPt);
364 fh3TrackPt->Fill(jetPt,invM,trackPt);
367 if(jetPt>0) z = trackPt / jetPt;
369 if(z>0) xi = TMath::Log(1/z);
371 fh3Xi->Fill(jetPt,invM,xi);
372 fh3Z->Fill(jetPt,invM,z);
375 //___________________________________________________________________________________
376 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AddToOutput(TList* list) const
378 // add histos to list
382 list->Add(fh3TrackPt);
390 //_______________________________________________________________________________________________________
391 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const char* name,
392 Int_t nPt, Float_t ptMin, Float_t ptMax,
393 Int_t nPhi, Float_t phiMin, Float_t phiMax,
394 Int_t nInvMass, Float_t invMassMin, Float_t invMassMax)
402 ,fNBinsInvMass(nInvMass)
403 ,fInvMassMin(invMassMin)
404 ,fInvMassMax(invMassMax)
408 // default constructor
411 //____________________________________________________________________________________________________________________________________
412 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const AliFragFuncHistosPhiCorrInvMass& copy)
414 ,fNBinsPt(copy.fNBinsPt)
417 ,fNBinsPhi(copy.fNBinsPhi)
418 ,fPhiMin(copy.fPhiMin)
419 ,fPhiMax(copy.fPhiMax)
420 ,fNBinsInvMass(copy.fNBinsInvMass)
421 ,fInvMassMin(copy.fInvMassMin)
422 ,fInvMassMax(copy.fInvMassMax)
423 ,fh3PhiCorr(copy.fh3PhiCorr)
424 ,fNamePhiCorr(copy.fNamePhiCorr)
429 //________________________________________________________________________________________________________________________________________________________________________
430 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& o)
435 TObject::operator=(o);
436 fNBinsPt = o.fNBinsPt;
439 fNBinsPhi = o.fNBinsPhi;
442 fNBinsInvMass = o.fNBinsInvMass;
443 fInvMassMin = o.fInvMassMin;
444 fInvMassMax = o.fInvMassMax;
446 fh3PhiCorr = o.fh3PhiCorr;
447 fNamePhiCorr = o.fNamePhiCorr;
453 //_________________________________________________________________________________________
454 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::~AliFragFuncHistosPhiCorrInvMass()
458 if(fh3PhiCorr) delete fh3PhiCorr;
461 //__________________________________________________________________________
462 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::DefineHistos()
464 // book jet QA histos
466 fh3PhiCorr = new TH3F(Form("fh3PhiCorrIM%s", fNamePhiCorr.Data()),
467 Form("%s: p_{t} - #phi - m_{inv} distribution",fNamePhiCorr.Data()),
468 fNBinsPt, fPtMin, fPtMax,
469 fNBinsPhi, fPhiMin, fPhiMax,
470 fNBinsInvMass, fInvMassMin, fInvMassMax);
472 AliAnalysisTaskJetChem::SetProperties(fh3PhiCorr, "p_{t} (GeV/c)", "#phi", "m_{inv} (GeV/c^2)");
475 //___________________________________________________________________________________________________________
476 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::FillPhiCorr(Float_t pt, Float_t phi, Float_t invM)
478 // fill jet QA histos
480 fh3PhiCorr->Fill(pt, phi, invM);
483 //______________________________________________________________________________________________
484 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AddToOutput(TList* list) const
486 // add histos to list
488 list->Add(fh3PhiCorr);
491 //____________________________________________________
492 void AliAnalysisTaskJetChem::UserCreateOutputObjects()
494 // create output objects
496 if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()");
498 // create list of tracks and jets
500 fTracksRecCuts = new TList();
501 fTracksRecCuts->SetOwner(kFALSE);
503 fJetsRecCuts = new TList();
504 fJetsRecCuts->SetOwner(kFALSE);
506 fListK0s = new TList();
507 fListK0s->SetOwner(kFALSE);
510 // Create histograms / output container
514 fCommonHistList = new TList();
516 Bool_t oldStatus = TH1::AddDirectoryStatus();
517 TH1::AddDirectory(kFALSE);
520 // histograms inherited from AliAnalysisTaskFragmentationFunction
522 fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
523 fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
524 fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
525 fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
526 fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
527 fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
528 fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
530 fh1EvtCent = new TH1F("fh1EvtCent","centrality",100,0.,100.);
531 fh1VertexNContributors = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
532 fh1VertexZ = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
533 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
534 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
535 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
536 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
537 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
538 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
540 fh1nRecJetsCuts = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
543 // histograms jetChem proper
545 fh1EvtMult = new TH1F("fh1EvtMult","multiplicity",1200,0.,12000.);
546 fh1K0Mult = new TH1F("fh1K0Mult","K0 multiplicity",500,0.,500.);
547 fh1dPhiJetK0 = new TH1F("fh1dPhiJetK0","",640,-1,5.4);
551 fFFHistosRecCuts = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
552 fFFNBinsPt, fFFPtMin, fFFPtMax,
553 fFFNBinsXi, fFFXiMin, fFFXiMax,
554 fFFNBinsZ , fFFZMin , fFFZMax);
556 fV0QAK0 = new AliFragFuncQATrackHistos("V0QAK0",fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
557 fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
558 fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
559 fQATrackHighPtThreshold);
561 fFFHistosRecCutsK0Evt = new AliFragFuncHistos("RecCutsK0Evt", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
562 fFFNBinsPt, fFFPtMin, fFFPtMax,
563 fFFNBinsXi, fFFXiMin, fFFXiMax,
564 fFFNBinsZ , fFFZMin , fFFZMax);
567 fFFHistosIMK0AllEvt = new AliFragFuncHistosInvMass("K0AllEvt", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
568 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
569 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
570 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
571 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
573 fFFHistosIMK0Jet = new AliFragFuncHistosInvMass("K0Jet", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
574 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
575 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
576 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
577 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
580 fFFHistosIMK0Cone = new AliFragFuncHistosInvMass("K0Cone", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
581 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
582 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
583 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
584 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
586 fFFHistosPhiCorrIMK0 = new AliFragFuncHistosPhiCorrInvMass("K0",fPhiCorrIMNBinsPt, fPhiCorrIMPtMin, fPhiCorrIMPtMax,
587 fPhiCorrIMNBinsPhi, fPhiCorrIMPhiMin, fPhiCorrIMPhiMax,
588 fPhiCorrIMNBinsInvM , fPhiCorrIMInvMMin , fPhiCorrIMInvMMax);
591 fV0QAK0->DefineHistos();
592 fFFHistosRecCuts->DefineHistos();
593 fFFHistosRecCutsK0Evt->DefineHistos();
594 fFFHistosIMK0AllEvt->DefineHistos();
595 fFFHistosIMK0Jet->DefineHistos();
596 fFFHistosIMK0Cone->DefineHistos();
597 fFFHistosPhiCorrIMK0->DefineHistos();
599 const Int_t saveLevel = 5;
602 fCommonHistList->Add(fh1EvtSelection);
603 fCommonHistList->Add(fh1EvtCent);
604 fCommonHistList->Add(fh1VertexNContributors);
605 fCommonHistList->Add(fh1VertexZ);
606 fCommonHistList->Add(fh1Xsec);
607 fCommonHistList->Add(fh1Trials);
608 fCommonHistList->Add(fh1PtHard);
609 fCommonHistList->Add(fh1PtHardTrials);
610 fCommonHistList->Add(fh1nRecJetsCuts);
611 fCommonHistList->Add(fh1EvtMult);
612 fCommonHistList->Add(fh1K0Mult);
613 fCommonHistList->Add(fh1dPhiJetK0);
615 fV0QAK0->AddToOutput(fCommonHistList);
616 fFFHistosRecCuts->AddToOutput(fCommonHistList);
617 fFFHistosRecCutsK0Evt->AddToOutput(fCommonHistList);
619 fFFHistosIMK0AllEvt->AddToOutput(fCommonHistList);
620 fFFHistosIMK0Jet->AddToOutput(fCommonHistList);
621 fFFHistosIMK0Cone->AddToOutput(fCommonHistList);
622 fFFHistosPhiCorrIMK0->AddToOutput(fCommonHistList);
625 // =========== Switch on Sumw2 for all histos ===========
626 for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
627 TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
630 THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
631 if(hnSparse) hnSparse->Sumw2();
635 TH1::AddDirectory(oldStatus);
638 //_______________________________________________
639 void AliAnalysisTaskJetChem::UserExec(Option_t *)
642 // Called for each event
644 if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
646 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
647 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
648 if(inputHandler->IsEventSelected() & AliVEvent::kMB){
649 if(fDebug > 1) Printf(" Trigger Selection: event ACCEPTED ... ");
650 fh1EvtSelection->Fill(1.);
652 fh1EvtSelection->Fill(0.);
653 if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
654 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
655 PostData(1, fCommonHistList);
660 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
662 if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
665 fMCEvent = MCEvent();
667 if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
670 // get AOD event from input/ouput
671 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
672 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
673 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
674 if(fUseAODInputJets) fAODJets = fAOD;
675 if (fDebug > 1) Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
678 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
679 if( handler && handler->InheritsFrom("AliAODHandler") ) {
680 fAOD = ((AliAODHandler*)handler)->GetAOD();
682 if (fDebug > 1) Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
686 if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
687 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
688 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
689 fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
690 if (fDebug > 1) Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
694 if(fNonStdFile.Length()!=0){
695 // case we have an AOD extension - fetch the jets from the extended output
697 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
698 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
700 if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
705 Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
709 Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
714 // event selection *****************************************
716 Double_t centPercent = -1;
719 if(handler && handler->InheritsFrom("AliAODInputHandler")){
720 // since it is not supported by the helper task define own classes
721 centPercent = fAOD->GetHeader()->GetCentrality();
723 if(centPercent>10) cl = 2;
724 if(centPercent>30) cl = 3;
725 if(centPercent>50) cl = 4;
728 cl = AliAnalysisHelperJetTasks::EventClass();
729 if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); // OB added
733 // event not in selected event class, reject event
734 if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
735 fh1EvtSelection->Fill(2.);
736 PostData(1, fCommonHistList);
741 // *** vertex cut ***
742 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
743 Int_t nTracksPrim = primVtx->GetNContributors();
744 fh1VertexNContributors->Fill(nTracksPrim);
746 if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
748 if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
749 fh1EvtSelection->Fill(3.);
750 PostData(1, fCommonHistList);
754 fh1VertexZ->Fill(primVtx->GetZ());
756 if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
757 if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
758 fh1EvtSelection->Fill(4.);
759 PostData(1, fCommonHistList);
763 TString primVtxName(primVtx->GetName());
765 if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
766 if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
767 fh1EvtSelection->Fill(5.);
768 PostData(1, fCommonHistList);
772 if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
773 fh1EvtSelection->Fill(0.);
774 fh1EvtCent->Fill(centPercent);
777 //___ get MC information __________________________________________________________________
779 Double_t ptHard = 0.;
780 Double_t nTrials = 1; // trials for MC trigger weight for real data
783 AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
784 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
785 AliGenHijingEventHeader* hijingGenHeader = 0x0;
788 if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
789 nTrials = pythiaGenHeader->Trials();
790 ptHard = pythiaGenHeader->GetPtHard();
792 fh1PtHard->Fill(ptHard);
793 fh1PtHardTrials->Fill(ptHard,nTrials);
796 } else { // no pythia, hijing?
798 if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
800 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
801 if(!hijingGenHeader){
802 Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
804 if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
808 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
812 //____ fetch jets ______________________________________________________________
814 Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);
815 Int_t nRecJetsCuts = 0;
816 if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
817 if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
818 if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
819 fh1nRecJetsCuts->Fill(nRecJetsCuts);
822 //____ fetch particles __________________________________________________________
824 Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);
825 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
826 if(fTracksRecCuts->GetEntries() != nTCuts)
827 Printf("%s:%d Mismatch selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
828 fh1EvtMult->Fill(fTracksRecCuts->GetEntries());
830 Int_t nK0s = GetListOfK0s(fListK0s,fK0Type);
831 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
832 if(nK0s != fListK0s->GetEntries()) Printf("%s:%d Mismatch selected K0s: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
833 fh1K0Mult->Fill(fListK0s->GetEntries());
835 // ___ V0 QA + K0 pt spectra all events _______________________________________________
837 if(fListK0s->GetEntries()>0){
839 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
841 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
844 Float_t trackPt = v0->Pt();
845 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
846 Double_t invM = v0->MassK0Short();
847 Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
849 fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt());
850 fFFHistosIMK0AllEvt->FillFF(trackPt, invM, jetPt, incrementJetPt);
855 //____ fill FF histos __________________________________________________________
857 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
859 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
861 if(ij==0){ // leading jet
863 TList* jettracklist = new TList();
866 if(GetFFRadius()<=0){
867 GetJetTracksTrackrefs(jettracklist, jet);
869 GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt);
872 for(Int_t it=0; it<jettracklist->GetSize(); ++it){
874 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
875 if(!trackVP)continue;
877 Float_t jetPt = jet->Pt();
878 Float_t trackPt = trackVP->Pt();
880 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
882 fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);
883 if(nK0s>0) fFFHistosRecCutsK0Evt->FillFF(trackPt, jetPt, incrementJetPt);
891 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // jet loop
893 AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(ij));
896 if(ij==0){ // leading jet
898 //fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
900 Float_t jetPt = jet->Pt();
902 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
904 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
909 TVector3 v0MomVect(v0Mom);
911 Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
913 Float_t trackPt = v0->Pt();
914 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
915 Double_t invM = v0->MassK0Short();
917 fFFHistosIMK0Jet->FillFF(trackPt, invM, jetPt, incrementJetPt);
918 fFFHistosPhiCorrIMK0->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetK0),invM);
920 if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
921 fh1dPhiJetK0->Fill(dPhiJetK0);
924 if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum
926 Bool_t incrementJetPt = kTRUE;
927 fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
931 TList* jetConeK0list = new TList();
934 GetJetTracksPointing(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPt);
936 if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
938 for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop K0s in jet cone
940 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
943 Double_t invM = v0->MassK0Short();
944 Float_t trackPt = v0->Pt();
945 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
947 //std::cout<<" trackPt "<<trackPt<<" invM "<<invM<<std::endl;
948 fFFHistosIMK0Cone->FillFF(trackPt, invM, jetPt, incrementJetPt);
950 if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum
952 Bool_t incrementJetPt = kTRUE;
953 fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
956 delete jetConeK0list;
960 fTracksRecCuts->Clear();
961 fJetsRecCuts->Clear();
965 PostData(1, fCommonHistList);
968 // ____________________________________________________________________________________________
969 void AliAnalysisTaskJetChem::SetProperties(TH3F* h,const char* x, const char* y, const char* z)
971 //Set properties of histos (x,y and z title)
976 h->GetXaxis()->SetTitleColor(1);
977 h->GetYaxis()->SetTitleColor(1);
978 h->GetZaxis()->SetTitleColor(1);
981 // ____________________________________________________________________________________________________________________________________________
982 Bool_t AliAnalysisTaskJetChem::IsAccepteddEdx(const Double_t mom,const Double_t signal, AliPID::EParticleType n, const Double_t cutnSig) const{
984 // apply TPC dE/dx cut similar as in AliTPCpidESD
985 // note: AliTPCidESD uses inner track param for momentum - not avaiable on AOD,
986 // so we use global track momentum
987 // should use separate parametrisation for MC and data, but probably ALEPH param & 7% resolution used here anyway not the last word
990 const Double_t kBBMIP(50.);
991 const Double_t kBBRes(0.07);
992 //const Double_t kBBRange(5.);
993 const Double_t kBBp1(0.76176e-1);
994 const Double_t kBBp2(10.632);
995 const Double_t kBBp3(0.13279e-4);
996 const Double_t kBBp4(1.8631);
997 const Double_t kBBp5(1.9479);
999 Double_t mass=AliPID::ParticleMass(n);
1000 Double_t betaGamma = mom/mass;
1002 const Float_t kmeanCorrection =0.1;
1003 Double_t bb = AliExternalTrackParam::BetheBlochAleph(betaGamma,kBBp1,kBBp2,kBBp3,kBBp4,kBBp5);
1004 Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
1005 Double_t bethe = bb * meanCorrection; // expected
1006 Double_t sigma = bethe * kBBRes;
1009 Double_t dedx = signal/kBBMIP; // measured
1011 Double_t nSig = (TMath::Abs(dedx - bethe))/sigma;
1013 if(nSig > cutnSig) return kFALSE;
1018 //___________________________________________________________________
1019 Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
1021 // K0 mass ? Use FF histo limits
1023 if(fFFIMInvMMin <= mass && mass <= fFFIMInvMMax) return kTRUE;
1028 //_____________________________________________________________________________________
1029 Int_t AliAnalysisTaskJetChem::GetListOfK0s(TList *list, const Int_t type)
1031 // fill list of V0s selected according to type
1034 if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
1038 if(type==kTrackUndef) return 0;
1040 for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
1042 AliAODv0* v0 = fAOD->GetV0(i);
1044 Bool_t isOnFly = v0->GetOnFlyStatus();
1046 if(!isOnFly && (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue;
1047 if( isOnFly && (type == kOffl || type == kOfflPID || type == kOffldEdx || type == kOfflPrim)) continue;
1049 Double_t massK0 = v0->MassK0Short();
1051 if(!(IsK0InvMass(massK0))) continue; // moved invMass cut for HI - otherwise too slow
1053 if(type == kOnFlyPID || type == kOfflPID){
1055 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1056 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1058 // AOD pid - cuts strongly into signal
1060 AliAODTrack::AODTrkPID_t mpPIDNeg = trackNeg->GetMostProbablePID();
1061 AliAODTrack::AODTrkPID_t mpPIDPos = trackPos->GetMostProbablePID();
1063 if(!( (mpPIDNeg == AliAODTrack::kPion) && (mpPIDPos == AliAODTrack::kPion) ) ) continue;
1066 if(type == kOnFlydEdx || type == kOffldEdx){
1068 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1069 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1071 AliAODPid* aodPidPos = trackPos->GetDetPid();
1072 AliAODPid* aodPidNeg = trackNeg->GetDetPid();
1074 Double_t dEdxPos = aodPidPos->GetTPCsignal();
1075 Double_t dEdxNeg = aodPidNeg->GetTPCsignal();
1077 Double_t momPos = trackPos->P();
1078 Double_t momNeg = trackNeg->P();
1080 Int_t cutnSigdEdx = 2;
1081 if(! (IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,cutnSigdEdx)) && (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,cutnSigdEdx)) ) continue;
1085 if(type == kOnFlyPrim || type == kOfflPrim){
1087 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1088 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1090 //std::cout<<" filer map trackPos "<<trackPos->GetFilterMap()<<" trackNeg "<<trackNeg->GetFilterMap()<<std::endl;
1092 if((fFilterMaskK0>0) && !(trackPos->TestFilterBit(fFilterMaskK0))) continue;
1093 if((fFilterMaskK0>0) && !(trackNeg->TestFilterBit(fFilterMaskK0))) continue;
1099 Int_t nK0s = list->GetSize();