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()
231 if(fListK0s) delete fListK0s;
234 //________________________________________________________________________________________________________________________________
235 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const char* name,
236 Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,
237 Int_t nInvMass, Float_t invMassMin, Float_t invMassMax,
238 Int_t nPt, Float_t ptMin, Float_t ptMax,
239 Int_t nXi, Float_t xiMin, Float_t xiMax,
240 Int_t nZ , Float_t zMin , Float_t zMax )
245 ,fNBinsInvMass(nInvMass)
246 ,fInvMassMin(invMassMin)
247 ,fInvMassMax(invMassMax)
263 // default constructor
267 //______________________________________________________________________________________________________________
268 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const AliFragFuncHistosInvMass& copy)
270 ,fNBinsJetPt(copy.fNBinsJetPt)
271 ,fJetPtMin(copy.fJetPtMin)
272 ,fJetPtMax(copy.fJetPtMax)
273 ,fNBinsInvMass(copy.fNBinsInvMass)
274 ,fInvMassMin(copy.fInvMassMin)
275 ,fInvMassMax(copy.fInvMassMax)
276 ,fNBinsPt(copy.fNBinsPt)
279 ,fNBinsXi(copy.fNBinsXi)
282 ,fNBinsZ(copy.fNBinsZ)
285 ,fh3TrackPt(copy.fh3TrackPt)
288 ,fh1JetPt(copy.fh1JetPt)
289 ,fNameFF(copy.fNameFF)
294 //______________________________________________________________________________________________________________________________________________________________________
295 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& o)
300 TObject::operator=(o);
301 fNBinsJetPt = o.fNBinsJetPt;
302 fJetPtMin = o.fJetPtMin;
303 fJetPtMax = o.fJetPtMax;
304 fNBinsInvMass = o.fNBinsInvMass;
305 fInvMassMin = o.fInvMassMin;
306 fInvMassMax = o.fInvMassMax;
307 fNBinsPt = o.fNBinsPt;
310 fNBinsXi = o.fNBinsXi;
316 fh3TrackPt = o.fh3TrackPt;
319 fh1JetPt = o.fh1JetPt;
326 //___________________________________________________________________________
327 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::~AliFragFuncHistosInvMass()
331 if(fh1JetPt) delete fh1JetPt;
332 if(fh3TrackPt) delete fh3TrackPt;
333 if(fh3Xi) delete fh3Xi;
334 if(fh3Z) delete fh3Z;
337 //_________________________________________________________________
338 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::DefineHistos()
342 fh1JetPt = new TH1F(Form("fh1FFJetPtIM%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
343 fh3TrackPt = new TH3F(Form("fh3FFTrackPtIM%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsPt, fPtMin, fPtMax);
344 fh3Xi = new TH3F(Form("fh3FFXiIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsXi, fXiMin, fXiMax);
345 fh3Z = new TH3F(Form("fh3FFZIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsZ, fZMin, fZMax);
347 AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{t} (GeV/c)", "entries");
348 AliAnalysisTaskJetChem::SetProperties(fh3TrackPt,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","p_{t} (GeV/c)");
349 AliAnalysisTaskJetChem::SetProperties(fh3Xi,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","#xi");
350 AliAnalysisTaskJetChem::SetProperties(fh3Z,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","z");
353 //________________________________________________________________________________________________________________________________
354 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::FillFF(Float_t trackPt, Float_t invM, Float_t jetPt, Bool_t incrementJetPt)
358 if(incrementJetPt) fh1JetPt->Fill(jetPt);
359 fh3TrackPt->Fill(jetPt,invM,trackPt);
362 if(jetPt>0) z = trackPt / jetPt;
364 if(z>0) xi = TMath::Log(1/z);
366 fh3Xi->Fill(jetPt,invM,xi);
367 fh3Z->Fill(jetPt,invM,z);
370 //___________________________________________________________________________________
371 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AddToOutput(TList* list) const
373 // add histos to list
377 list->Add(fh3TrackPt);
385 //_______________________________________________________________________________________________________
386 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const char* name,
387 Int_t nPt, Float_t ptMin, Float_t ptMax,
388 Int_t nPhi, Float_t phiMin, Float_t phiMax,
389 Int_t nInvMass, Float_t invMassMin, Float_t invMassMax)
397 ,fNBinsInvMass(nInvMass)
398 ,fInvMassMin(invMassMin)
399 ,fInvMassMax(invMassMax)
403 // default constructor
406 //____________________________________________________________________________________________________________________________________
407 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const AliFragFuncHistosPhiCorrInvMass& copy)
409 ,fNBinsPt(copy.fNBinsPt)
412 ,fNBinsPhi(copy.fNBinsPhi)
413 ,fPhiMin(copy.fPhiMin)
414 ,fPhiMax(copy.fPhiMax)
415 ,fNBinsInvMass(copy.fNBinsInvMass)
416 ,fInvMassMin(copy.fInvMassMin)
417 ,fInvMassMax(copy.fInvMassMax)
418 ,fh3PhiCorr(copy.fh3PhiCorr)
419 ,fNamePhiCorr(copy.fNamePhiCorr)
424 //________________________________________________________________________________________________________________________________________________________________________
425 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& o)
430 TObject::operator=(o);
431 fNBinsPt = o.fNBinsPt;
434 fNBinsPhi = o.fNBinsPhi;
437 fNBinsInvMass = o.fNBinsInvMass;
438 fInvMassMin = o.fInvMassMin;
439 fInvMassMax = o.fInvMassMax;
441 fh3PhiCorr = o.fh3PhiCorr;
442 fNamePhiCorr = o.fNamePhiCorr;
448 //_________________________________________________________________________________________
449 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::~AliFragFuncHistosPhiCorrInvMass()
453 if(fh3PhiCorr) delete fh3PhiCorr;
456 //__________________________________________________________________________
457 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::DefineHistos()
459 // book jet QA histos
461 fh3PhiCorr = new TH3F(Form("fh3PhiCorrIM%s", fNamePhiCorr.Data()),
462 Form("%s: p_{t} - #phi - m_{inv} distribution",fNamePhiCorr.Data()),
463 fNBinsPt, fPtMin, fPtMax,
464 fNBinsPhi, fPhiMin, fPhiMax,
465 fNBinsInvMass, fInvMassMin, fInvMassMax);
467 AliAnalysisTaskJetChem::SetProperties(fh3PhiCorr, "p_{t} (GeV/c)", "#phi", "m_{inv} (GeV/c^2)");
470 //___________________________________________________________________________________________________________
471 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::FillPhiCorr(Float_t pt, Float_t phi, Float_t invM)
473 // fill jet QA histos
475 fh3PhiCorr->Fill(pt, phi, invM);
478 //______________________________________________________________________________________________
479 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AddToOutput(TList* list) const
481 // add histos to list
483 list->Add(fh3PhiCorr);
486 //____________________________________________________
487 void AliAnalysisTaskJetChem::UserCreateOutputObjects()
489 // create output objects
491 if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()");
493 // create list of tracks and jets
495 fTracksRecCuts = new TList();
496 fTracksRecCuts->SetOwner(kFALSE);
498 fJetsRecCuts = new TList();
499 fJetsRecCuts->SetOwner(kFALSE);
501 fListK0s = new TList();
502 fListK0s->SetOwner(kFALSE);
505 // Create histograms / output container
509 fCommonHistList = new TList();
511 Bool_t oldStatus = TH1::AddDirectoryStatus();
512 TH1::AddDirectory(kFALSE);
515 // histograms inherited from AliAnalysisTaskFragmentationFunction
517 fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
518 fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
519 fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
520 fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
521 fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
522 fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
523 fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
525 fh1EvtCent = new TH1F("fh1EvtCent","centrality",100,0.,100.);
526 fh1VertexNContributors = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
527 fh1VertexZ = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
528 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
529 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
530 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
531 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
532 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
533 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
535 fh1nRecJetsCuts = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
538 // histograms jetChem proper
540 fh1EvtMult = new TH1F("fh1EvtMult","multiplicity",1200,0.,12000.);
541 fh1K0Mult = new TH1F("fh1K0Mult","K0 multiplicity",500,0.,500.);
542 fh1dPhiJetK0 = new TH1F("fh1dPhiJetK0","",640,-1,5.4);
546 fFFHistosRecCuts = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
547 fFFNBinsPt, fFFPtMin, fFFPtMax,
548 fFFNBinsXi, fFFXiMin, fFFXiMax,
549 fFFNBinsZ , fFFZMin , fFFZMax);
551 fV0QAK0 = new AliFragFuncQATrackHistos("V0QAK0",fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
552 fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
553 fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
554 fQATrackHighPtThreshold);
556 fFFHistosRecCutsK0Evt = new AliFragFuncHistos("RecCutsK0Evt", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
557 fFFNBinsPt, fFFPtMin, fFFPtMax,
558 fFFNBinsXi, fFFXiMin, fFFXiMax,
559 fFFNBinsZ , fFFZMin , fFFZMax);
562 fFFHistosIMK0AllEvt = new AliFragFuncHistosInvMass("K0AllEvt", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
563 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
564 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
565 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
566 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
568 fFFHistosIMK0Jet = new AliFragFuncHistosInvMass("K0Jet", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
569 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
570 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
571 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
572 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
575 fFFHistosIMK0Cone = new AliFragFuncHistosInvMass("K0Cone", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
576 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
577 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
578 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
579 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
581 fFFHistosPhiCorrIMK0 = new AliFragFuncHistosPhiCorrInvMass("K0",fPhiCorrIMNBinsPt, fPhiCorrIMPtMin, fPhiCorrIMPtMax,
582 fPhiCorrIMNBinsPhi, fPhiCorrIMPhiMin, fPhiCorrIMPhiMax,
583 fPhiCorrIMNBinsInvM , fPhiCorrIMInvMMin , fPhiCorrIMInvMMax);
586 fV0QAK0->DefineHistos();
587 fFFHistosRecCuts->DefineHistos();
588 fFFHistosRecCutsK0Evt->DefineHistos();
589 fFFHistosIMK0AllEvt->DefineHistos();
590 fFFHistosIMK0Jet->DefineHistos();
591 fFFHistosIMK0Cone->DefineHistos();
592 fFFHistosPhiCorrIMK0->DefineHistos();
594 const Int_t saveLevel = 5;
597 fCommonHistList->Add(fh1EvtSelection);
598 fCommonHistList->Add(fh1EvtCent);
599 fCommonHistList->Add(fh1VertexNContributors);
600 fCommonHistList->Add(fh1VertexZ);
601 fCommonHistList->Add(fh1Xsec);
602 fCommonHistList->Add(fh1Trials);
603 fCommonHistList->Add(fh1PtHard);
604 fCommonHistList->Add(fh1PtHardTrials);
605 fCommonHistList->Add(fh1nRecJetsCuts);
606 fCommonHistList->Add(fh1EvtMult);
607 fCommonHistList->Add(fh1K0Mult);
608 fCommonHistList->Add(fh1dPhiJetK0);
610 fV0QAK0->AddToOutput(fCommonHistList);
611 fFFHistosRecCuts->AddToOutput(fCommonHistList);
612 fFFHistosRecCutsK0Evt->AddToOutput(fCommonHistList);
614 fFFHistosIMK0AllEvt->AddToOutput(fCommonHistList);
615 fFFHistosIMK0Jet->AddToOutput(fCommonHistList);
616 fFFHistosIMK0Cone->AddToOutput(fCommonHistList);
617 fFFHistosPhiCorrIMK0->AddToOutput(fCommonHistList);
620 // =========== Switch on Sumw2 for all histos ===========
621 for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
622 TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
625 THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
626 if(hnSparse) hnSparse->Sumw2();
630 TH1::AddDirectory(oldStatus);
633 //_______________________________________________
634 void AliAnalysisTaskJetChem::UserExec(Option_t *)
637 // Called for each event
639 if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
641 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
642 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
643 if(inputHandler->IsEventSelected() & AliVEvent::kMB){
644 if(fDebug > 1) Printf(" Trigger Selection: event ACCEPTED ... ");
645 fh1EvtSelection->Fill(1.);
647 fh1EvtSelection->Fill(0.);
648 if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
649 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
650 PostData(1, fCommonHistList);
655 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
657 if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
660 fMCEvent = MCEvent();
662 if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
665 // get AOD event from input/ouput
666 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
667 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
668 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
669 if(fUseAODInputJets) fAODJets = fAOD;
670 if (fDebug > 1) Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
673 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
674 if( handler && handler->InheritsFrom("AliAODHandler") ) {
675 fAOD = ((AliAODHandler*)handler)->GetAOD();
677 if (fDebug > 1) Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
681 if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
682 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
683 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
684 fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
685 if (fDebug > 1) Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
689 if(fNonStdFile.Length()!=0){
690 // case we have an AOD extension - fetch the jets from the extended output
692 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
693 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
695 if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
700 Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
704 Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
709 // event selection *****************************************
711 Double_t centPercent = -1;
714 if(handler && handler->InheritsFrom("AliAODInputHandler")){
715 // since it is not supported by the helper task define own classes
716 centPercent = fAOD->GetHeader()->GetCentrality();
718 if(centPercent>10) cl = 2;
719 if(centPercent>30) cl = 3;
720 if(centPercent>50) cl = 4;
723 cl = AliAnalysisHelperJetTasks::EventClass();
724 if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); // OB added
728 // event not in selected event class, reject event
729 if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
730 fh1EvtSelection->Fill(2.);
731 PostData(1, fCommonHistList);
736 // *** vertex cut ***
737 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
738 Int_t nTracksPrim = primVtx->GetNContributors();
739 fh1VertexNContributors->Fill(nTracksPrim);
741 if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
743 if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
744 fh1EvtSelection->Fill(3.);
745 PostData(1, fCommonHistList);
749 fh1VertexZ->Fill(primVtx->GetZ());
751 if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
752 if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
753 fh1EvtSelection->Fill(4.);
754 PostData(1, fCommonHistList);
758 TString primVtxName(primVtx->GetName());
760 if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
761 if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
762 fh1EvtSelection->Fill(5.);
763 PostData(1, fCommonHistList);
767 if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
768 fh1EvtSelection->Fill(0.);
769 fh1EvtCent->Fill(centPercent);
772 //___ get MC information __________________________________________________________________
774 Double_t ptHard = 0.;
775 Double_t nTrials = 1; // trials for MC trigger weight for real data
778 AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
779 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
780 AliGenHijingEventHeader* hijingGenHeader = 0x0;
783 if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
784 nTrials = pythiaGenHeader->Trials();
785 ptHard = pythiaGenHeader->GetPtHard();
787 fh1PtHard->Fill(ptHard);
788 fh1PtHardTrials->Fill(ptHard,nTrials);
791 } else { // no pythia, hijing?
793 if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
795 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
796 if(!hijingGenHeader){
797 Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
799 if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
803 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
807 //____ fetch jets ______________________________________________________________
809 Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);
810 Int_t nRecJetsCuts = 0;
811 if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
812 if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
813 if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
814 fh1nRecJetsCuts->Fill(nRecJetsCuts);
817 //____ fetch particles __________________________________________________________
819 Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);
820 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
821 if(fTracksRecCuts->GetEntries() != nTCuts)
822 Printf("%s:%d Mismatch selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
823 fh1EvtMult->Fill(fTracksRecCuts->GetEntries());
825 Int_t nK0s = GetListOfK0s(fListK0s,fK0Type);
826 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
827 if(nK0s != fListK0s->GetEntries()) Printf("%s:%d Mismatch selected K0s: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
828 fh1K0Mult->Fill(fListK0s->GetEntries());
830 // ___ V0 QA + K0 pt spectra all events _______________________________________________
832 if(fListK0s->GetEntries()>0){
834 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
836 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
839 Float_t trackPt = v0->Pt();
840 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
841 Double_t invM = v0->MassK0Short();
842 Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
844 fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt());
845 fFFHistosIMK0AllEvt->FillFF(trackPt, invM, jetPt, incrementJetPt);
850 //____ fill FF histos __________________________________________________________
852 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
854 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
855 Double_t jetPt = jet->Pt();
857 if(ij==0){ // leading jet
859 TList* jettracklist = new TList();
861 Bool_t isBadJet = kFALSE;
863 if(GetFFRadius()<=0){
864 GetJetTracksTrackrefs(jettracklist, jet, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
866 GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
869 if(GetFFMinNTracks()>0 && jettracklist->GetSize() <= GetFFMinNTracks()) isBadJet = kTRUE;
870 if(isBadJet) continue;
873 for(Int_t it=0; it<jettracklist->GetSize(); ++it){
875 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
876 if(!trackVP)continue;
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);
890 // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
892 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
894 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
899 TVector3 v0MomVect(v0Mom);
901 Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
903 Float_t trackPt = v0->Pt();
904 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
905 Double_t invM = v0->MassK0Short();
907 fFFHistosIMK0Jet->FillFF(trackPt, invM, jetPt, incrementJetPt);
908 fFFHistosPhiCorrIMK0->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetK0),invM);
910 if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
911 fh1dPhiJetK0->Fill(dPhiJetK0);
915 if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum
917 Bool_t incrementJetPt = kTRUE;
918 fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
922 TList* jetConeK0list = new TList();
923 Double_t sumPtK0 = 0.;
924 Bool_t isBadJetK0 = kFALSE; // dummy, do not use
926 GetJetTracksPointing(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPtK0, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetK0);
929 if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
931 for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop K0s in jet cone
933 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
936 Double_t invM = v0->MassK0Short();
937 Float_t trackPt = v0->Pt();
938 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
940 //std::cout<<" trackPt "<<trackPt<<" invM "<<invM<<std::endl;
941 fFFHistosIMK0Cone->FillFF(trackPt, invM, jetPt, incrementJetPt);
943 if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum
945 Bool_t incrementJetPt = kTRUE;
946 fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
949 delete jetConeK0list;
953 fTracksRecCuts->Clear();
954 fJetsRecCuts->Clear();
958 PostData(1, fCommonHistList);
961 // ____________________________________________________________________________________________
962 void AliAnalysisTaskJetChem::SetProperties(TH3F* h,const char* x, const char* y, const char* z)
964 //Set properties of histos (x,y and z title)
969 h->GetXaxis()->SetTitleColor(1);
970 h->GetYaxis()->SetTitleColor(1);
971 h->GetZaxis()->SetTitleColor(1);
974 // ____________________________________________________________________________________________________________________________________________
975 Bool_t AliAnalysisTaskJetChem::IsAccepteddEdx(const Double_t mom,const Double_t signal, AliPID::EParticleType n, const Double_t cutnSig) const{
977 // apply TPC dE/dx cut similar as in AliTPCpidESD
978 // note: AliTPCidESD uses inner track param for momentum - not avaiable on AOD,
979 // so we use global track momentum
980 // should use separate parametrisation for MC and data, but probably ALEPH param & 7% resolution used here anyway not the last word
983 const Double_t kBBMIP(50.);
984 const Double_t kBBRes(0.07);
985 //const Double_t kBBRange(5.);
986 const Double_t kBBp1(0.76176e-1);
987 const Double_t kBBp2(10.632);
988 const Double_t kBBp3(0.13279e-4);
989 const Double_t kBBp4(1.8631);
990 const Double_t kBBp5(1.9479);
992 Double_t mass=AliPID::ParticleMass(n);
993 Double_t betaGamma = mom/mass;
995 const Float_t kmeanCorrection =0.1;
996 Double_t bb = AliExternalTrackParam::BetheBlochAleph(betaGamma,kBBp1,kBBp2,kBBp3,kBBp4,kBBp5);
997 Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
998 Double_t bethe = bb * meanCorrection; // expected
999 Double_t sigma = bethe * kBBRes;
1002 Double_t dedx = signal/kBBMIP; // measured
1004 Double_t nSig = (TMath::Abs(dedx - bethe))/sigma;
1006 if(nSig > cutnSig) return kFALSE;
1011 //___________________________________________________________________
1012 Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
1014 // K0 mass ? Use FF histo limits
1016 if(fFFIMInvMMin <= mass && mass <= fFFIMInvMMax) return kTRUE;
1021 //_____________________________________________________________________________________
1022 Int_t AliAnalysisTaskJetChem::GetListOfK0s(TList *list, const Int_t type)
1024 // fill list of V0s selected according to type
1027 if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
1031 if(type==kTrackUndef) return 0;
1033 for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
1035 AliAODv0* v0 = fAOD->GetV0(i);
1037 Bool_t isOnFly = v0->GetOnFlyStatus();
1039 if(!isOnFly && (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue;
1040 if( isOnFly && (type == kOffl || type == kOfflPID || type == kOffldEdx || type == kOfflPrim)) continue;
1042 Double_t massK0 = v0->MassK0Short();
1044 if(!(IsK0InvMass(massK0))) continue; // moved invMass cut for HI - otherwise too slow
1046 if(type == kOnFlyPID || type == kOfflPID){
1048 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1049 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1051 // AOD pid - cuts strongly into signal
1053 AliAODTrack::AODTrkPID_t mpPIDNeg = trackNeg->GetMostProbablePID();
1054 AliAODTrack::AODTrkPID_t mpPIDPos = trackPos->GetMostProbablePID();
1056 if(!( (mpPIDNeg == AliAODTrack::kPion) && (mpPIDPos == AliAODTrack::kPion) ) ) continue;
1059 if(type == kOnFlydEdx || type == kOffldEdx){
1061 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1062 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1064 AliAODPid* aodPidPos = trackPos->GetDetPid();
1065 AliAODPid* aodPidNeg = trackNeg->GetDetPid();
1067 Double_t dEdxPos = aodPidPos->GetTPCsignal();
1068 Double_t dEdxNeg = aodPidNeg->GetTPCsignal();
1070 Double_t momPos = trackPos->P();
1071 Double_t momNeg = trackNeg->P();
1073 Int_t cutnSigdEdx = 2;
1074 if(! (IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,cutnSigdEdx)) && (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,cutnSigdEdx)) ) continue;
1078 if(type == kOnFlyPrim || type == kOfflPrim){
1080 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1081 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1083 //std::cout<<" filer map trackPos "<<trackPos->GetFilterMap()<<" trackNeg "<<trackNeg->GetFilterMap()<<std::endl;
1085 if((fFilterMaskK0>0) && !(trackPos->TestFilterBit(fFilterMaskK0))) continue;
1086 if((fFilterMaskK0>0) && !(trackNeg->TestFilterBit(fFilterMaskK0))) continue;
1092 Int_t nK0s = list->GetSize();