]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetChem.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding jet chemistry, not yet included in .pkg
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetChem.cxx
CommitLineData
9b2de807 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id:$ */
17
18
19#include <TChain.h>
20#include <TFile.h>
21#include <TList.h>
22#include <TH1F.h>
23#include <TH2F.h>
24#include <TProfile.h>
25#include <TVector3.h>
26#include <TLorentzVector.h>
27#include <TMath.h>
28#include <TTree.h>
29#include <TDatabasePDG.h>
30
31#include "AliAnalysisTaskJetChem.h"
32#include "AliAnalysisManager.h"
33#include "AliMCEventHandler.h"
34#include "AliAODInputHandler.h"
35#include "AliAODHandler.h"
36#include "AliAODJet.h"
37#include "AliAODTrack.h"
38#include "AliAODMCParticle.h"
39
40#include "AliGenPythiaEventHeader.h"
41#include "AliAnalysisHelperJetTasks.h"
42#include "AliESDtrack.h"
43#include "AliExternalTrackParam.h"
44
45//
46// Analysis class for jet chemistry studies
47// based on AliAnalysisTaskUE by Arian Abrahantes Quintana and Enesto Lopez
48// contact: Oliver Busch, o.busch@gsi.de
49
50// This class needs as input AOD with track and Jets
51// the output is a list of histograms
52//
53// AOD can be either connected to the InputEventHandler
54// for a chain of AOD files
55// or
56// to the OutputEventHandler
57// for a chain of ESD files, so this case class should be
58// in the train after the Jet finder
59//
60
61
62
63
64ClassImp(AliAnalysisTaskJetChem)
65
66////////////////////////////////////////////////////////////////////////
67
68
69//____________________________________________________________________
70AliAnalysisTaskJetChem:: AliAnalysisTaskJetChem(const char* name): AliAnalysisTaskSE(name),
71fDebug(kFALSE),
72fDeltaAOD(kFALSE),
73fDeltaAODBranch(""),
74fAODBranch("jets"),
75fDeltaAODBranchMC(""),
76fAODBranchMC("jetsMC"),
77fArrayJetsAOD(0x0),
78fArrayJetsMC(0x0),
79fAOD(0x0),
80fAODjets(0x0),
81fListOfHistos(0x0),
82fJetsOnFly(kFALSE),
83fUseLOConeJets(kFALSE),
84fUseLOConeMCJets(kFALSE),
85fUsePythiaJets(kFALSE),
86fConeRadius(0.),
87fTrackPtCutJF(0),
88fFilterBitJF(0x01),
89fRequireITSRefitJF(kFALSE),
90fRejectK0TracksJF(kFALSE),
91fJetPtCut(5.0),
92fJetEtaCut(0.5),
93fFilterBit(0x10),
94fTrackPtCut(0.),
95fTrackEtaCut(0.9),
96fUseOnFlyV0s(kFALSE),
97fCutnSigdEdx(0.),
98fUseAODMCTracksForUE(kFALSE),
99fAreaReg(1.0),
100fAvgTrials(1),
101fhPrimVertexRho(0x0),
102fhPrimVertexZ(0x0),
103fhNJets(0x0),
104fhNJetsMC(0x0),
105fhLeadingEta(0x0),
106fhLeadingPhi(0x0),
107fhLeadingPt(0x0),
108fhLeadingPtDiffr(0x0),
109fhLeadingEtaMC(0x0),
110fhLeadingPhiMC(0x0),
111fhLeadingPtMC(0x0),
112fhLeadingPtMCDiffr(0x0),
113fhEtaTracks(0x0),
114fhPhiTracks(0x0),
115fhPtTracks(0x0),
116fhEtaMCTracks(0x0),
117fhPhiMCTracks(0x0),
118fhPtMCTracks(0x0),
119fhnTracksVsPtLeading(0x0),
120fhdNdEtaPhiDist(0x0),
121fhRegionSumPtMaxVsEt(0x0),
122fhRegionMultMaxVsEt(0x0),
123fhRegionSumPtMinVsEt(0x0),
124fhRegionMultMinVsEt(0x0),
125fhNV0s(0x0),
126fhV0onFly(0x0),
127fhV0DCADaughters(0x0),
128fhV0Radius(0x0),
129fhV0DCAToVertex(0x0),
130fhV0DCAToVertexK0(0x0),
131fhV0InvMassK0(0x0),
132fhV0InvMassK0JetEvt(0x0),
133fhV0InvMassLambda(0x0),
134fhV0InvMassAntiLambda(0x0),
135fhV0InvMassLambdaJetEvt(0x0),
136fhV0InvMassAntiLambdaJetEvt(0x0),
137fhdROpanK0VsPt(0x0),
138fhdPhiJetV0(0x0),
139fhdPhiJetK0(0x0),
140fhdRJetK0(0x0),
141fhdNdptV0(0x0),
142fhdNdptK0(0x0),
143fhPtVsEtaK0(0x0),
144fhV0InvMassK0DCA(0x0),
145fhV0InvMassK0DCAdEdx(0x0),
146fhV0InvMassK0DCAPID(0x0),
147fhV0InvMassLambdaDCAdEdx(0x0),
148fhV0InvMassAntiLambdaDCAdEdx(0x0),
149fhdNdptK0DCA(0x0),
150fhdNdptK0DCAdEdx(0x0),
151fhV0InvMassK0Min(0x0),
152fhV0InvMassLambdaMin(0x0),
153fhV0InvMassAntiLambdaMin(0x0),
154fhV0InvMassK0Max(0x0),
155fhV0InvMassLambdaMax(0x0),
156fhV0InvMassAntiLambdaMax(0x0),
157fhV0InvMassK0Jet(0x0),
158fhV0InvMassLambdaJet(0x0),
159fhV0InvMassAntiLambdaJet(0x0),
160fhV0InvMassK0Lambda(0x0),
161fhdNdptK0JetEvt(0x0),
162fhdNdzK0(0x0),
163fhdNdxiK0(0x0),
164fhdNdzK0Max(0x0),
165fhdNdxiK0Max(0x0),
166fhdNdzLambdaMax(0x0),
167fhdNdxiLambdaMax(0x0),
168fhdNdptK0Max(0x0),
169fhdNdptLambdaMax(0x0),
170fhdNdzK0Min(0x0),
171fhdNdxiK0Min(0x0),
172fhdNdzLambdaMin(0x0),
173fhdNdxiLambdaMin(0x0),
174fhdNdptK0Min(0x0),
175fhdNdptLambdaMin(0x0),
176fhdNdzK0Jet(0x0),
177fhdNdxiK0Jet(0x0),
178fhdNdzLambdaJet(0x0),
179fhdNdxiLambdaJet(0x0),
180fhdNdptK0Jet(0x0),
181fhdNdptLambdaJet(0x0),
182fhdEdxVsMomV0(0x0),
183fhdEdxVsMomV0pidEdx(0x0),
184fhdEdxVsMomV0piPID(0x0),
185fhdPhiJetK0MC(0x0),
186fhdRJetK0MC(0x0),
187fhdNdptK0MCMax(0x0),
188fhdNdptLambdaLambdaBarMCMax(0x0),
189fhdNdptK0MCMin(0x0),
190fhdNdptLambdaLambdaBarMCMin(0x0),
191fhdNdptOmegaMCMin(0x0),
192fhdNdptOmegaBarMCMin(0x0),
193fhdNdptchPiMCMin(0x0),
194fhdNdptchKMCMin(0x0),
195fhdNdptpMCMin(0x0),
196fhdNdptpBarMCMin(0x0),
197fhdNdptK0MCJet(0x0),
198fhdNdptLambdaLambdaBarMCJet(0x0),
199fhdNdptchPiMCJet(0x0),
200fhdNdptchKMCJet(0x0),
201fhdNdptpMCJet(0x0),
202fhdNdptpBarMCJet(0x0),
203fhdRV0MC(0x0),
204fhPIDMC(0x0),
205fhPIDMCAll(0x0),
206fhPIDMCMin(0x0),
207fhPIDMCJet(0x0),
208fhdNdptgammaMC(0x0),
209fhdNdptchPiMC(0x0),
210fhdNdptpi0MC(0x0),
211fhdNdptK0MC(0x0),
212fhdNdptchKMC(0x0),
213fhdNdptpMC(0x0),
214fhdNdptpBarMC(0x0),
215fhdNdptLambdaLambdaBarMC(0x0),
216fhdNdptOmegaMC(0x0),
217fhdNdptOmegaBarMC(0x0),
218fhdNdxiMC(0x0),
219fhdNdxiK0MC(0x0),
220fhdNdxiK0MCJet(0x0),
221fhdNdzK0MC(0x0),
222fhdNdzK0MCJet(0x0),
223fhdNdptK0MCJetEvt(0x0),
224fhnJetsAODvsMC(0x0),
225fhLeadingPtAODvsMC(0x0),
226fhLeadingEtaAODvsMC(0x0),
227fhLeadingPhiAODvsMC(0x0),
228fhnTracksLeadingAODvsMC(0x0),
229fhLeadingdRAODMC(0x0),
230fhLeadingPtAODvsMCdRcut(0x0),
231fhdnTracksVsdPtLeadingAODMC(0x0),
232fhnTracksJetVsPtAOD(0x0),
233fhnTracksJetVsPtAODquarkEv(0x0),
234fhRadiusJetVsPtAOD(0x0),
235fhnTracksJetVsPtMC(0x0),
236fhnTracksJetVsPtMCquarkEv(0x0),
237fhRadiusJetVsPtMC(0x0),
238fhnTracksJetVsPtMCK0(0x0),
239fhnTracksJetVsPtMCK0quarkEv(0x0),
240fhRadiusJetVsPtMCK0(0x0),
241fhnTracksJetVsPtAODK0(0x0),
242fhnTracksJetVsPtAODK0quarkEv(0x0),
243fhRadiusJetVsPtAODK0(0x0),
244fhnTracksJetVsPtAODpKch(0x0),
245fhRadiusJetVsPtAODpKch(0x0),
246fhPythiaProcess(0x0),
247fhPythiaProcessK0(0x0),
248fhPythiaProcessKch(0x0),
249fhPythiaProcessp(0x0),
250fhPythiaProcesspbar(0x0),
251fh1Xsec(0x0),
252fh1Trials(0x0)
253{
254 // Default constructor
255
256 fAreaReg = 2*TMath::Pi()/6.0 * 2*fTrackEtaCut;
257 fpdgdb = TDatabasePDG::Instance();
258
259
260 // Output slot #1 writes into a TList container, 0 reserved for std AOD output
261 DefineOutput(1, TList::Class());
262}
263
264//______________________________________________________________
265Bool_t AliAnalysisTaskJetChem::UserNotify()
266{
267 //
268 // read the cross sections
269 // and number of trials from pyxsec.root
270 //
271
272 fAvgTrials = 1;
273 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
274 Float_t xsection = 0;
275 Float_t ftrials = 1;
276 if(tree){
277 TFile *curfile = tree->GetCurrentFile();
278 if (!curfile) {
279 Error("Notify","No current file");
280 return kFALSE;
281 }
282 if(!fh1Xsec||!fh1Trials){
283 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
284 return kFALSE;
285 }
286 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
287 fh1Xsec->Fill("<#sigma>",xsection);
288 // construct average trials
289 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
290 if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries;
291 }
292 return kTRUE;
293}
294
295//____________________________________________________________________
296void AliAnalysisTaskJetChem::UserCreateOutputObjects()
297{
298 // Create the output container
299 //
300 AliInfo("UserCreateOutPutObjects()");
301 //
302 // Histograms
303
304 OpenFile(1);
305 CreateHistos();
306 PostData(1, fListOfHistos); // PostData at least once
307
308}
309
310//____________________________________________________________________
311void AliAnalysisTaskJetChem::UserExec(Option_t */*option*/)
312{
313 // Execute analysis for current event
314 //
315 if (fDebug > 3) AliInfo( " Processing event..." );
316
317 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
318
319 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
320 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
321 if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler");
322 // Case when jets are reconstructed on the fly from AOD tracks
323 // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
324 // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
325 // different parameters to default ones stored in the AOD or to use a different algorithm
326 if( fJetsOnFly ) {
327 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
328 if( handler && handler->InheritsFrom("AliAODHandler") ) {
329 fAODjets = ((AliAODHandler*)handler)->GetAOD();
330 if (fDebug > 1) AliInfo(" ==== Jets from AliAODHandler");
331 }
332 } else {
333 fAODjets = fAOD;
334 if (fDebug > 1) AliInfo(" ==== Jets from AliAODInputHandler");
335 }
336 } else {
337 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
338 if( handler && handler->InheritsFrom("AliAODHandler") ) {
339 fAOD = ((AliAODHandler*)handler)->GetAOD();
340 fAODjets = fAOD;
341 if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
342 } else {
343 AliFatal("I can't get any AOD Event Handler");
344 return;
345 }
346 }
347
348 // -------------
349
350 // fetch the pythia header info and get the trials
351 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
352 Float_t nTrials = 1;
353 if (mcHandler) {
354 AliMCEvent* mcEvent = mcHandler->MCEvent();
355 if (mcEvent) {
356 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
357 if(pythiaGenHeader){
358 nTrials = pythiaGenHeader->Trials();
359 }
360 }
361 }
362
363 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
364
365 AnalyseEvent();
366
367 // Post the data
368 PostData(1, fListOfHistos);
369}
370
371//____________________________________________________________________
372void AliAnalysisTaskJetChem::AnalyseEvent()
373{
374 // check Trigger
375
376 TString firedTriggerClasses = fAOD->GetFiredTriggerClasses();
377 // AliInfo(Form("firedTriggerClasses %s",firedTriggerClasses.Data()));
378 // if(firedTriggerClasses.Length() > 0 && !firedTriggerClasses.Contains("CINT1B")) return;
379
380 AliAODVertex *primVertex = fAOD->GetPrimaryVertex(); // is this from tracks or from SPD ? SPD should only be used if no track vertex
381 if(!primVertex){
382 AliInfo(Form("no primVertex found - skip event"));
383 return;
384 }
385
386 Double_t vertX = primVertex->GetX();
387 Double_t vertY = primVertex->GetY();
388 Double_t vertZ = primVertex->GetZ();
389
390 Double_t vertRho = TMath::Sqrt(vertX*vertX+vertY*vertY);
391
392 Int_t nVertContributors = primVertex->GetNContributors();
393 if(nVertContributors <1){
394 AliInfo(Form("prim vertex no contributors - skip event"));
395 return;
396 }
397
398 fhPrimVertexRho->Fill(vertRho);
399 fhPrimVertexZ->Fill(vertZ);
400
401 Int_t pythiaPID = GetPythiaProcessID();
402
403 // ------------------------------------------------
404 // Find Leading Jets 1,2,3
405 // (could be skipped if Jets are sort by Pt...)
406
407 //Int_t index1 = -1;
408 //Int_t nTracksLeading = 0;
409
410 Int_t nJetsAOD = 0;
411 AliAODJet* leadingJetAOD = NULL; // non-zero if leading jet survives acc, pt cut
412 Int_t indexLeadingAOD = -1;
413 Double_t ptLeadingAOD = 0.;
414 Int_t indexMaxRegionAOD = 0; // initialize with '0', '+-1' transverse regions
415
416 Int_t nJetsMC = 0;
417 AliAODJet* leadingJetMC = NULL; // non-zero if leading jet survives acc, pt cut
418 Int_t indexLeadingMC = -1;
419 Double_t ptLeadingMC = 0.;
420 Int_t indexMaxRegionMC = 0; // initialize with '0', '+-1' transverse regions
421
422
423 if(fUseLOConeJets) fArrayJetsAOD = FindChargedParticleJets();
424 else{ // Use jets on AOD/DeltaAOD
425
426 if(fDeltaAOD){
427 if (fDebug > 1) AliInfo(" ==== Jets From Delta-AODs !");
428 if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fDeltaAODBranch.Data()));
429 fArrayJetsAOD = (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranch.Data());
430 if (!fArrayJetsAOD){
431 AliFatal(" No jet-array! ");
432 return;
433 }
434 }
435 else{
436 if (fDebug > 1) AliInfo(" ==== AOD jets: Read Standard-AODs !");
437 if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fAODBranch.Data()));
438
439 fArrayJetsAOD = (TClonesArray*)fAODjets->FindListObject(fAODBranch.Data());
440 }
441 }
442
443
444 if(fUseLOConeMCJets) fArrayJetsMC = FindChargedParticleJetsMC();
445 else if(fUsePythiaJets) fArrayJetsMC = GetPythiaJets();
446 else{
447
448 if(fDeltaAOD){
449 if (fDebug > 1) AliInfo(" ==== MC Jets From Delta-AODs !");
450 if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fDeltaAODBranchMC.Data()));
451 fArrayJetsMC = (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranchMC.Data());
452 if (!fArrayJetsMC){
453 AliFatal(" No jet-array! ");
454 return;
455 }
456 }
457 else{
458 if (fDebug > 1) AliInfo(" ==== MC jets: Read Standard-AODs !");
459 if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fAODBranchMC.Data()));
460
461 fArrayJetsMC = (TClonesArray*)fAODjets->FindListObject(fAODBranchMC.Data());
462 }
463 }
464
465 if(fArrayJetsAOD) nJetsAOD = fArrayJetsAOD->GetEntries();
466 if(fArrayJetsMC) nJetsMC = fArrayJetsMC->GetEntries();
467
468 fhNJets->Fill(nJetsAOD);
469 fhNJetsMC->Fill(nJetsMC);
470 fhnJetsAODvsMC->Fill(nJetsMC,nJetsAOD);
471
472 if(fDebug>1) AliInfo(Form("AOD %d jets",nJetsAOD));
473
474
475 // find leading jet AOD
476 for(Int_t i=0; i<nJetsAOD; ++i){
477
478 AliAODJet* jet = (AliAODJet*) fArrayJetsAOD->At(i);
479 Double_t jetPt = jet->Pt();
480 Double_t jetEta = jet->Eta();
481
482 if((jetPt > ptLeadingAOD) && (TMath::Abs(jetEta)<fJetEtaCut)){
483 ptLeadingAOD = jetPt;
484 indexLeadingAOD = i;
485 }
486 }
487
488 // find leading jet MC
489 for(Int_t i=0; i<nJetsMC; ++i){
490
491 AliAODJet* jet = (AliAODJet*) fArrayJetsMC->At(i);
492 Double_t jetPt = jet->Pt();
493 Double_t jetEta = jet->Eta();
494
495 if((jetPt > ptLeadingMC) && (TMath::Abs(jetEta)<fJetEtaCut)){
496 ptLeadingMC = jetPt;
497 indexLeadingMC = i;
498 }
499 }
500
501
502 //cout<<" ptLeadingAOD "<<ptLeadingAOD<<" MC "<<ptLeadingMC<<endl;
503
504 if(indexLeadingAOD>=0){ // event with jet
505
506 Double_t etaLeadingAOD = ((AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD))->Eta();
507 Double_t phiLeadingAOD = ((AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD))->Phi();
508 Int_t nTracksLeadingAOD = ((AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD))->GetRefTracks()->GetEntriesFast();
509
510 if(fDebug>1) AliInfo(Form("\n Pt Leading AOD Jet = %6.1f eta=%5.3f, nTracks %d",ptLeadingAOD,etaLeadingAOD,nTracksLeadingAOD));
511
512 fhLeadingEta->Fill(etaLeadingAOD);
513
514 if(TMath::Abs(etaLeadingAOD)<fJetEtaCut){ // leading jet eta cut
515
516 fhnTracksVsPtLeading->Fill(ptLeadingAOD,nTracksLeadingAOD);
517 fhLeadingPt->Fill(ptLeadingAOD);
518 if(IsDiffractiveEvent(pythiaPID)) fhLeadingPtDiffr->Fill(ptLeadingAOD);
519
520 if(ptLeadingAOD>fJetPtCut){ // leading jet pt cut
521
522 leadingJetAOD = (AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD);
523
524 fhLeadingPhi->Fill(TMath::ATan2(TMath::Sin(phiLeadingAOD),TMath::Cos(phiLeadingAOD))); // -pi to pi
525
526 // ----------------------------------------------
527 // Find max and min regions
528 Double_t sumPtRegionPosit = 0.;
529 Double_t sumPtRegionNegat = 0.;
530 Int_t nTrackRegionPosit = 0;
531 Int_t nTrackRegionNegat = 0;
532
533 Int_t nTracks = fAOD->GetNTracks();
534
535 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
536
537 AliAODTrack* part = fAOD->GetTrack(ipart);
538 if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
539 if (!part->IsPrimaryCandidate()) continue; // reject whatever is not linked to collision point
540 // PID Selection: Reject everything but hadrons
541 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
542 part->GetMostProbablePID()==AliAODTrack::kKaon ||
543 part->GetMostProbablePID()==AliAODTrack::kProton;
544 if (!isHadron ) continue;
545 if (!part->Charge() ) continue; //Only charged
546 if (part->Pt() < fTrackPtCut ) continue;
547 if(TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
548
549 TVector3 partVect(part->Px(), part->Py(), part->Pz());
550
551 Int_t region = IsTrackInsideRegion(leadingJetAOD,&partVect );
552
553 if (region > 0) {
554 sumPtRegionPosit += part->Pt();
555 nTrackRegionPosit++;
556 }
557 if (region < 0) {
558 sumPtRegionNegat += part->Pt();
559 nTrackRegionNegat++;
560 }
561 } // tracks loop
562
563 // fill sumPt and mult density for transverse regions
564
565 if( sumPtRegionPosit > sumPtRegionNegat ) {
566 FillSumPtRegion( ptLeadingAOD, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
567 }
568 else {
569 FillSumPtRegion( ptLeadingAOD, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
570 }
571 if (nTrackRegionPosit > nTrackRegionNegat ) {
572 FillMultRegion( ptLeadingAOD, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg);
573 }
574 else {
575 FillMultRegion( ptLeadingAOD, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg);
576 }
577
578 indexMaxRegionAOD = (sumPtRegionPosit > sumPtRegionNegat) ? 1 : -1;
579
580 } // leading jet pt cut
581 } // leading jet eta cut
582 } // jet event
583
584
585 if(indexLeadingMC>=0){ // event with MC jet
586
587 Double_t etaLeadingMC = ((AliAODJet*) fArrayJetsMC->At(indexLeadingMC))->Eta();
588 Double_t phiLeadingMC = ((AliAODJet*) fArrayJetsMC->At(indexLeadingMC))->Phi();
589 Int_t nTracksLeadingMC = ((AliAODJet*) fArrayJetsMC->At(indexLeadingMC))->GetRefTracks()->GetEntriesFast();
590
591 if(fDebug>1) AliInfo(Form("\n Pt Leading MC Jet = %6.1f eta=%5.3f, nTracks %d",ptLeadingMC,etaLeadingMC,nTracksLeadingMC));
592
593 fhLeadingEtaMC->Fill(etaLeadingMC);
594
595 if(TMath::Abs(etaLeadingMC)<fJetEtaCut){ // leading jet eta cut
596
597 fhLeadingPtMC->Fill(ptLeadingMC);
598 if(IsDiffractiveEvent(pythiaPID)) fhLeadingPtMCDiffr->Fill(ptLeadingMC);
599
600 if(ptLeadingMC>fJetPtCut){ // leading jet pt cut
601
602 leadingJetMC = (AliAODJet*) fArrayJetsMC->At(indexLeadingMC);
603
604 fhLeadingPhiMC->Fill(TMath::ATan2(TMath::Sin(phiLeadingMC),TMath::Cos(phiLeadingMC))); // -pi to pi
605
606 // ----------------------------------------------
607 // Find max and min regions
608 Double_t sumPtRegionPosit = 0;
609 Double_t sumPtRegionNegat = 0;
610 Int_t nTrackRegionPosit = 0;
611 Int_t nTrackRegionNegat = 0;
612
613 TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
614
615 Int_t ntrks = farray->GetEntries();
616 if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
617
618 for(Int_t i =0 ; i < ntrks; i++){
619
620 AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
621 //Cuts
622 if (!(mctrk->IsPhysicalPrimary())) continue;
623 //if (!(mctrk->IsPrimary())) continue;
624
625 if (mctrk->Charge() == 0 || mctrk->Charge()==-99) continue;
626
627 if (mctrk->Pt() < fTrackPtCut ) continue;
628 if( TMath::Abs(mctrk->Eta()) > fTrackEtaCut ) continue;
629
630 Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
631 TMath::Abs(mctrk->GetPdgCode())==2212 ||
632 TMath::Abs(mctrk->GetPdgCode())==321;
633
634 if (!isHadron) continue;
635
636 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
637
638 Int_t region = IsTrackInsideRegion(leadingJetMC,&partVect );
639
640 if (region > 0) {
641 sumPtRegionPosit += mctrk->Pt();
642 nTrackRegionPosit++;
643 }
644 if (region < 0) {
645 sumPtRegionNegat += mctrk->Pt();
646 nTrackRegionNegat++;
647 }
648 } // AliAODMCParticle loop
649
650 indexMaxRegionMC = (sumPtRegionPosit > sumPtRegionNegat) ? 1 : -1;
651
652 } // leading jet pt
653 } // leading jet eta
654 } // found jet
655
656
657 Bool_t foundK0AOD = kFALSE; // event with K0 ?
658 Bool_t foundK0MC = kFALSE; // event with K0 ?
659
660 CheckV0s(leadingJetAOD,indexMaxRegionAOD,foundK0AOD); // here leadingJetAOD/MC nonzero if jet passes eta & pt cut
661 CheckMCParticles(leadingJetMC,indexMaxRegionMC,foundK0MC);
662 CompLeadingJets(leadingJetAOD,leadingJetMC,pythiaPID,foundK0AOD,foundK0MC);
663
664 if(fUseLOConeJets && fArrayJetsAOD){
665 fArrayJetsAOD->Delete(); // no 'Clear': AliAODjet contains TMomentum and TRefArray
666 }
667 if(fUseLOConeMCJets && fArrayJetsMC){
668 fArrayJetsMC->Delete(); // no 'Clear': AliAODjet contains TMomentum and TRefArray
669 }
670}
671
672// __________________________________________________________________
673
674
675Double_t AliAnalysisTaskJetChem::GetJetRadius(const AliAODJet* jet, const Double_t energyFrac){
676
677 // calc jet radius containing fraction energyFrac of full jet pt
678
679 Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
680
681 Double_t deltaR[nTracks];
682 Double_t pt[nTracks];
683 Int_t index[nTracks];
684 for(int i=0; i<nTracks; i++) index[i] = -999;
685
686 Double_t ptTot = 0;
687
688 for(int i=0; i<nTracks; i++){
689
690 AliAODTrack* track = (AliAODTrack*) jet->GetRefTracks()->At(i);
691
692 TLorentzVector *mom4Jet = jet->MomentumVector();
693 TVector3 mom3Jet = mom4Jet->Vect();
694
695 Double_t trackMom[3];
696 track->PxPyPz(trackMom);
697 TVector3 mom3Track(trackMom);
698
699 Double_t dR = mom3Jet.DeltaR(mom3Track);
700
701 deltaR[i] = dR;
702 pt[i] = track->Pt();
703
704 ptTot += pt[i];
705 }
706
707 //cout<<" ptTot "<<ptTot<<" jetPt "<<jet->Pt()<<endl; // Xcheck
708
709 TMath::Sort(nTracks,deltaR,index,kFALSE); // sort in decreasing order
710
711 Double_t ptSum = 0;
712
713 for(int i=0; i<nTracks; i++){
714
715 Int_t ind = index[i];
716
717 Double_t ptTrack = pt[ind];
718 Double_t dR = deltaR[ind];
719
720 ptSum += ptTrack;
721
722 if(ptSum >= ptTot*energyFrac) return dR;
723 }
724
725 return -1;
726
727}
728
729//____________________________________________________________________
730void AliAnalysisTaskJetChem::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
731{
732 // Fill sumPt of control regions
733
734 fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
735 fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
736}
737
738//____________________________________________________________________
739void AliAnalysisTaskJetChem::FillMultRegion(Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin)
740{
741 // Fill Nch multiplicity of control regions
742
743 fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
744 fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
745
746}
747
748//____________________________________________________________________
749Int_t AliAnalysisTaskJetChem::IsTrackInsideRegion(const AliAODJet* aodjetVect,const TVector3 *partVect)
750{
751 // return de region in delta phi
752 // -1 negative delta phi
753 // 1 positive delta phi
754 // 0 outside region
755
756 TLorentzVector* jetVectLorentz = aodjetVect->MomentumVector();
757 TVector3 jetVect = jetVectLorentz->Vect();
758
759
760 static const Double_t k60rad = 60.*TMath::Pi()/180.;
761 static const Double_t k120rad = 120.*TMath::Pi()/180.;
762
763 Int_t region = 0;
764 if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
765 // transverse regions
766 if (jetVect.DeltaPhi(*partVect) < -k60rad && jetVect.DeltaPhi(*partVect) > -k120rad ) region = -1;
767 if (jetVect.DeltaPhi(*partVect) > k60rad && jetVect.DeltaPhi(*partVect) < k120rad ) region = 1;
768
769 return region;
770}
771
772//____________________________________________________________________
773TClonesArray* AliAnalysisTaskJetChem::FindChargedParticleJetsMC()
774{
775 // CDF jet finder:
776 // loop over pt-ordered list of tracks, combine tracks within jet cone, recalc cone axis after each step
777 // based on implementation by Arian Abrahantes Quintana and Enesto Lopez
778
779 TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
780 if(!farray){
781 AliInfo("no mcparticles branch");
782 return 0;
783 }
784
785 Int_t nTracks = farray->GetEntries();
786
787 if( !nTracks ) return 0;
788 TObjArray tracks(nTracks);
789
790 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
791
792 AliAODMCParticle* part = (AliAODMCParticle*)farray->At(ipart);
793
794 if(!part->IsPhysicalPrimary()) continue;
795
796 // exclude neutrinos
797 Int_t pdg = TMath::Abs(part->GetPdgCode());
798 if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE;
799
800 if( !part->Charge() ) continue;
801 fhEtaMCTracks->Fill(part->Eta());
802
803 if(TMath::Abs(part->Eta()) > fTrackEtaCut) continue;
804 fhPtMCTracks->Fill(part->Pt());
805
806 if( part->Pt() < fTrackPtCutJF ) continue;
807 fhPhiMCTracks->Fill(TMath::ATan2(TMath::Sin(part->Phi()),TMath::Cos(part->Phi())));
808
809 tracks.AddLast(part);
810
811 }
812
813 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
814
815 nTracks = tracks.GetEntriesFast();
816
817 if( !nTracks ) return 0;
818 TObjArray *jets = new TObjArray(nTracks);
819 TIter itrack(&tracks);
820 while( nTracks ) {
821 // 2- Start with the highest pT particle ...
822 Float_t px,py,pz,pt;
823
824 AliAODMCParticle* track = (AliAODMCParticle*)itrack.Next();
825
826 if( !track ) continue;
827 px = track->Px();
828 py = track->Py();
829 pz = track->Pz();
830 pt = track->Pt();
831 jets->AddLast( new AliAODJet(px,py,pz,pt) ); // Use the energy member to store Pt
832 AliAODJet* jet = (AliAODJet*)jets->Last();
833 jet->AddTrack(track);
834 tracks.Remove(track);
835
836
837 TVector2 jetVect2(jet->Px(),jet->Py());
838
839 // 3- Go to the next highest pT particle not already included...
840 AliAODMCParticle* track1;
841 while ( (track1 = (AliAODMCParticle*)(itrack.Next())) ) {
842 TVector2 trackVect2(track1->Px(),track1->Py());
843 Double_t dphi = trackVect2.DeltaPhi(jetVect2);
844
845 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
846 dphi*dphi );
847
848 if( r < fConeRadius ) {
849 Double_t fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
850 // recalculating the centroid
851 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
852
853 jetVect2.SetMagPhi(jet->E()/fPt,jetVect2.Phi());
854 trackVect2.SetMagPhi(track1->Pt()/fPt,trackVect2.Phi());
855
856 TVector2 sumVect2 = jetVect2+trackVect2;
857 Double_t phi = sumVect2.Phi();
858
859 //jet->SetPtEtaPhiE( 1., eta, phi, fPt );
860 ((TLorentzVector*) jet->MomentumVector())->SetPtEtaPhiE(1,eta,phi,fPt);
861
862 jet->AddTrack(track1);
863 tracks.Remove(track1);
864 }
865 }
866
867 tracks.Compress();
868
869 nTracks = tracks.GetEntries();
870 // 4- Continue until all particles are in a jet.
871 itrack.Reset();
872 } // end while nTracks
873
874 // Convert to AODjets....
875 Int_t njets = jets->GetEntriesFast();
876
877 TClonesArray* aodjets = new TClonesArray("AliAODJet",njets);
878 aodjets->SetOwner(kTRUE);
879
880 Int_t count = 0;
881 for(Int_t ijet=0; ijet<njets; ++ijet) {
882 AliAODJet* jet = (AliAODJet*)jets->At(ijet);
883
884 if (jet->E() < fJetPtCut) continue;
885 Float_t px, py,pz,en; // convert to 4-vector
886 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
887 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
888 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
889 en = TMath::Sqrt(px * px + py * py + pz * pz);
890 jet->SetPxPyPzE(px,py,pz,en);
891
892
893 TClonesArray &tmpaodjets = *aodjets;
894 new(tmpaodjets[count++]) AliAODJet(*jet);
895 //aodjets->AddLast( new AliAODJet(*jet));
896 //aodjets->AddLast( new AliAODJet(px, py, pz, en) );
897 }
898 // Order jets according to their pT .
899 QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
900
901 // debug
902 //if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
903
904 jets->Delete(); // OB - should I cleanup or leave it to garbage collection ?
905 delete jets;
906
907 return aodjets;
908}
909
910// ___________________________________________________________________
911
912Bool_t AliAnalysisTaskJetChem::IsTrackFromK0(const Int_t indexTrack){
913
914 // check wether track with index is from K0 (= V0 with proper inv mass)
915
916 TClonesArray* tracks = fAOD->GetTracks();
917
918 for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
919
920 AliAODv0* v0 = fAOD->GetV0(i);
921
922 Bool_t isOnFly = v0->GetOnFlyStatus();
923
924 if( (fUseOnFlyV0s && !isOnFly) || (!fUseOnFlyV0s && isOnFly) ) continue;
925
926 Double_t massK0 = v0->MassK0Short();
927
928 // FIXME: here should go cuts (at least V0 quality)
929
930 if(IsK0InvMass(massK0)){
931
932 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
933 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
934
935 Int_t indexPos = tracks->IndexOf(trackPos);
936 Int_t indexNeg = tracks->IndexOf(trackNeg);
937
938
939 if(indexPos == indexTrack){ return kTRUE;}
940 if(indexNeg == indexTrack){ return kTRUE;}
941
942 }
943 }
944
945 return kFALSE;
946
947}
948
949//____________________________________________________________________
950TClonesArray* AliAnalysisTaskJetChem::FindChargedParticleJets()
951{
952
953 // CDF jet finder:
954 // loop over pt-ordered list of tracks, combine tracks within jet cone, recalc cone axis after each step
955 // based on implementation by Arian Abrahantes Quintana and Enesto Lopez
956 // ref: PHYSICAL REVIEW D 65 092002, CDF Collaboration
957
958
959 // 1 - Order all charged particles according to their pT .
960 Int_t nTracks = fAOD->GetNTracks();
961
962 if( !nTracks ) return 0;
963 TObjArray tracks(nTracks);
964
965 // standard cuts + ITS refit = stdrd PWG4 cut (compose them for productions with old ESDFilter task)
966
967 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
968 AliAODTrack* part = fAOD->GetTrack( ipart );
969
970 UInt_t status = part->GetStatus();
971
972 if( !part->TestFilterBit(fFilterBitJF) ) continue; // track cut selection
973 if(fRequireITSRefitJF && ((status&AliESDtrack::kITSrefit)==0)) continue;
974
975 if(!part->Charge() ) continue;
976 fhEtaTracks->Fill(part->Eta());
977
978 if(TMath::Abs(part->Eta()) > fTrackEtaCut) continue;
979
980 if(fRejectK0TracksJF && IsTrackFromK0(ipart)) continue;
981
982 fhPtTracks->Fill(part->Pt());
983
984 if( part->Pt() < fTrackPtCutJF ) continue;
985
986 fhPhiTracks->Fill(TMath::ATan2(TMath::Sin(part->Phi()),TMath::Cos(part->Phi())));
987
988 tracks.AddLast(part);
989 }
990
991 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
992
993 nTracks = tracks.GetEntriesFast();
994
995 if( !nTracks ) return 0;
996 TObjArray *jets = new TObjArray(nTracks);
997 TIter itrack(&tracks);
998 while( nTracks ) {
999 // 2- Start with the highest pT particle ...
1000 Float_t px,py,pz,pt;
1001 AliAODTrack* track = (AliAODTrack*)itrack.Next();
1002 if( !track ) continue;
1003 px = track->Px();
1004 py = track->Py();
1005 pz = track->Pz();
1006 pt = track->Pt();
1007 //jets->AddLast( new TLorentzVector(px, py, pz, pt) ); // Use the energy member to store Pt
1008 jets->AddLast( new AliAODJet(px,py,pz,pt) ); // Use the energy member to store Pt
1009 //TLorentzVector* jet = (TLorentzVector*)jets->Last();
1010 AliAODJet* jet = (AliAODJet*)jets->Last();
1011 jet->AddTrack(track);
1012 tracks.Remove( track );
1013
1014 TVector2 jetVect2(jet->Px(),jet->Py()); // OB
1015
1016 // 3- Go to the next highest pT particle not already included...
1017 AliAODTrack* track1;
1018 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
1019 //Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-track1->Phi()); // OB remove
1020 TVector2 trackVect2(track1->Px(),track1->Py());
1021 Double_t dphi = trackVect2.DeltaPhi(jetVect2);
1022
1023 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
1024 dphi*dphi );
1025
1026 if( r < fConeRadius ) {
1027 Double_t fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
1028 // recalculating the centroid
1029 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
1030
1031 //Double_t phi = jet->Phi()*jet->E()/fPt + track1->Phi()*track1->Pt()/fPt; // OB - remove
1032 // OB - recalc phi via weighted 2vectors
1033 jetVect2.SetMagPhi(jet->E()/fPt,jetVect2.Phi());
1034 trackVect2.SetMagPhi(track1->Pt()/fPt,trackVect2.Phi());
1035
1036 TVector2 sumVect2 = jetVect2+trackVect2;
1037 Double_t phi = sumVect2.Phi();
1038
1039 //jet->SetPtEtaPhiE( 1., eta, phi, fPt );
1040 ((TLorentzVector*) jet->MomentumVector())->SetPtEtaPhiE(1,eta,phi,fPt);
1041
1042 jet->AddTrack(track1);
1043 tracks.Remove(track1);
1044 }
1045 }
1046
1047 tracks.Compress();
1048
1049 nTracks = tracks.GetEntries();
1050 // 4- Continue until all particles are in a jet.
1051 itrack.Reset();
1052 } // end while nTracks
1053
1054 // Convert to AODjets....
1055 Int_t njets = jets->GetEntriesFast();
1056
1057 TClonesArray* aodjets = new TClonesArray("AliAODJet",njets);
1058 aodjets->SetOwner(kTRUE);
1059
1060 Int_t count = 0;
1061
1062 for(Int_t ijet=0; ijet<njets; ++ijet) {
1063 //TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
1064 AliAODJet* jet = (AliAODJet*)jets->At(ijet);
1065
1066 if (jet->E() < fJetPtCut) continue;
1067 Float_t px, py,pz,en; // convert to 4-vector
1068 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
1069 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
1070 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
1071 en = TMath::Sqrt(px * px + py * py + pz * pz);
1072 jet->SetPxPyPzE(px,py,pz,en);
1073
1074 TClonesArray &tmpaodjets = *aodjets;
1075 new(tmpaodjets[count++]) AliAODJet(*jet);
1076 //aodjets->AddLast( new AliAODJet(*jet));
1077 //aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1078 }
1079 // Order jets according to their pT .
1080 //QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() ); // not for TClonesArray
1081
1082 if (fDebug>3) AliInfo(Form(" %d Charged jets found - after cuts %d \n",njets,count));
1083
1084 jets->Delete(); // OB: cleanup
1085 delete jets;
1086
1087
1088 return aodjets;
1089}
1090
1091//____________________________________________________________________
1092
1093TClonesArray* AliAnalysisTaskJetChem::GetPythiaJets(){
1094
1095 // return Pythia jets in jet acc
1096
1097 // note: present verion of AliMCEventHandler expects "AliAOD.root" (not "AliAODs.root"), otherwise galice.root will not be found in proper dir
1098 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1099 if (!mcHandler) {
1100 Printf("ERROR: Could not retrieve MC event handler");
1101 return 0;
1102 }
1103
1104 AliMCEvent* mcEvent = mcHandler->MCEvent();
1105 if (!mcEvent) {
1106 Printf("ERROR: Could not retrieve MC event");
1107 return 0;
1108 }
1109
1110
1111 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1112 if(!pythiaGenHeader){
1113 Printf("ERROR: Could not retrieve pythiaGenHeader");
1114 return 0;
1115 }
1116
1117 //Get Jets from MC header
1118 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
1119
1120 TClonesArray* aodjets = new TClonesArray("AliAODJet",nPythiaGenJets);
1121 aodjets->SetOwner(kTRUE);
1122
1123 int count = 0;
1124 for(int ip = 0; ip<nPythiaGenJets; ++ip){
1125 Float_t p[4];
1126 pythiaGenHeader->TriggerJet(ip,p);
1127 TVector3 tempVect(p[0],p[1],p[2]);
1128 if ( TMath::Abs(tempVect.Eta())>fJetEtaCut ) continue;
1129
1130 TClonesArray &tmpaodjets = *aodjets;
1131 new(tmpaodjets[count++]) AliAODJet(p[0],p[1],p[2],p[3]);
1132 //aodjets->AddLast(new AliAODJet(p[0],p[1],p[2],p[3]));
1133 }
1134
1135 return aodjets;
1136}
1137
1138
1139//____________________________________________________________________
1140void AliAnalysisTaskJetChem::QSortTracks(TObjArray &a, Int_t first, Int_t last)
1141{
1142 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1143
1144 static TObject *tmp;
1145 static int i; // "static" to save stack space
1146 int j;
1147
1148 while (last - first > 1) {
1149 i = first;
1150 j = last;
1151 for (;;) {
1152 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
1153 ;
1154 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
1155 ;
1156 if (i >= j)
1157 break;
1158
1159 tmp = a[i];
1160 a[i] = a[j];
1161 a[j] = tmp;
1162 }
1163 if (j == first) {
1164 ++first;
1165 continue;
1166 }
1167 tmp = a[first];
1168 a[first] = a[j];
1169 a[j] = tmp;
1170 if (j - first < last - (j + 1)) {
1171 QSortTracks(a, first, j);
1172 first = j + 1; // QSortTracks(j + 1, last);
1173 } else {
1174 QSortTracks(a, j + 1, last);
1175 last = j; // QSortTracks(first, j);
1176 }
1177 }
1178}
1179
1180//------------------------------------------------------------------
1181
1182TH1F* AliAnalysisTaskJetChem::CreatePIDhisto(const char* name){
1183
1184 // create histogram
1185
1186 TH1F* result = new TH1F(name,"",50,0,50);
1187 result->SetOption("E");
1188
1189 // bin equal Geant ID
1190
1191 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(1),"photon");
1192 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(2),"e+");
1193 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(3),"e-");
1194 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(4),"e-neutrino");
1195 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(5),"mu+");
1196 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(6),"mu-");
1197 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(7),"pi0");
1198 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(8),"pi+");
1199 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(9),"pi-");
1200 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(10),"K long");
1201 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(11),"K+");
1202 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(12),"K-");
1203 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(13),"n");
1204 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(14),"p");
1205 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(15),"anti-proton");
1206 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(16),"K short");
1207 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(17),"eta");
1208 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(18),"Lambda");
1209 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(19),"Sigma+");
1210 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(20),"Sigma0");
1211 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(21),"Sigma-");
1212 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(22),"Xi0");
1213 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(23),"Xi-");
1214 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(24),"Omega-");
1215 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(25),"anti-neutron");
1216 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(26),"anti-Lambda");
1217 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(27),"anti-Sigma-");
1218 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(28),"anti-Sigma0");
1219 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(29),"anti-Sigma+");
1220 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(30),"anti-Xi0");
1221 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(31),"anti-Xi+");
1222 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(32),"anti-Omega+");
1223 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(33),"tau+");
1224 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(34),"tau-");
1225 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(35),"D+");
1226 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(36),"D-");
1227 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(37),"D0");
1228 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(38),"anti-D0");
1229 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(39),"Ds+");
1230 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(40),"anti Ds-");
1231 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(41),"Lamba_c+");
1232 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(42),"W+");
1233 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(43),"W-");
1234 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(44),"Z0");
1235 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(45),"d");
1236 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(46),"t");
1237 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(47),"alpha");
1238 result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(48),"G_nu");
1239
1240
1241 result->GetXaxis()->LabelsOption("v"); // "u" ?
1242
1243 return result;
1244}
1245
1246//------------------------------------------------------------------
1247
1248TH1F* AliAnalysisTaskJetChem::CreatePythiaIDhisto(const char* name){
1249
1250 // create histogram
1251
1252 TH1F* result = new TH1F(name,"",12,0,12);
1253 result->SetOption("E");
1254
1255 result->GetXaxis()->SetBinLabel(kPythiaPIDP11Bin,"qq #rightarrow qq"); // ISUB = 11
1256 result->GetXaxis()->SetBinLabel(kPythiaPIDP12Bin,"q#bar{q} #rightarrow q#bar{q}"); // ISUB = 12
1257 result->GetXaxis()->SetBinLabel(kPythiaPIDP13Bin,"q#bar{q} #rightarrow gg"); // ISUB = 13
1258 result->GetXaxis()->SetBinLabel(kPythiaPIDP28Bin,"qg #rightarrow qg "); // ISUB = 28
1259 result->GetXaxis()->SetBinLabel(kPythiaPIDP53Bin,"gg #rightarrow q#bar{q}"); // ISUB = 53
1260 result->GetXaxis()->SetBinLabel(kPythiaPIDP68Bin,"gg #rightarrow gg"); // ISUB = 68
1261
1262 result->GetXaxis()->LabelsOption("v"); // "u" ?
1263
1264 return result;
1265}
1266
1267
1268//------------------------------------------------------------------
1269
1270void AliAnalysisTaskJetChem::FillPythiaIDhisto(TH1F* h, const Int_t PID){
1271
1272 // fill Pyhtia PID histogram
1273
1274 Int_t bin = -1;
1275
1276 if(PID == 11) bin = kPythiaPIDP11Bin;
1277 else if(PID == 12) bin = kPythiaPIDP12Bin;
1278 else if(PID == 13) bin = kPythiaPIDP13Bin;
1279 else if(PID == 28) bin = kPythiaPIDP28Bin;
1280 else if(PID == 53) bin = kPythiaPIDP53Bin;
1281 else if(PID == 68) bin = kPythiaPIDP68Bin;
1282 else{
1283 if(PID != -1) AliInfo(Form("unknown PID %d",PID));
1284 return;
1285 }
1286
1287 h->Fill(h->GetBinCenter(bin),1);
1288}
1289
1290
1291//____________________________________________________________________
1292
1293
1294void AliAnalysisTaskJetChem::FillPIDhisto(TH1F* hist, Int_t pdg, Float_t weight){
1295
1296 // convert pdg code to Geant ID and fill corresponding bin
1297
1298 Int_t fGID = fpdgdb->ConvertPdgToGeant3(pdg);
1299
1300 //cout<<" pdg "<<pdg<<" GID "<<GID<<endl;
1301
1302 hist->Fill(fGID,weight);
1303
1304 if(fGID>hist->GetBinLowEdge(hist->GetNbinsX()+1) ||
1305 fGID<hist->GetBinLowEdge(1)){
1306
1307 AliError(Form("fGID %d for pdg %d exceeding histo limits ",fGID,pdg));
1308 }
1309
1310 if(fGID == 0){
1311 //AliError(Form("fGID 0 for pdg %d ",pdg));
1312 }
1313
1314}
1315
1316//____________________________________________________________________
1317void AliAnalysisTaskJetChem::CreateHistos(){
1318
1319 // make histos
1320
1321 fListOfHistos = new TList();
1322 fListOfHistos->SetOwner(kTRUE);
1323
1324 fhPrimVertexRho = new TH1F("hPrimVertexRho","",100,0,1);
1325 fhPrimVertexRho->SetXTitle("");
1326 fhPrimVertexRho->Sumw2();
1327 fListOfHistos->Add( fhPrimVertexRho );
1328
1329 fhPrimVertexZ = new TH1F("hPrimVertexZ","",200,0,20);
1330 fhPrimVertexZ->SetXTitle("");
1331 fhPrimVertexZ->Sumw2();
1332 fListOfHistos->Add( fhPrimVertexZ );
1333
1334 fhNJets = new TH1F("hNJets","",10, 0, 10);
1335 fhNJets->SetXTitle("# of jets");
1336 fhNJets->Sumw2();
1337 fListOfHistos->Add( fhNJets );
1338
1339 fhNJetsMC = new TH1F("hNJetsMC","",10,0,10);
1340 fhNJetsMC->SetXTitle("# of jets");
1341 fhNJetsMC->Sumw2();
1342 fListOfHistos->Add( fhNJetsMC );
1343
1344 fhLeadingEta = new TH1F("hLeadingEta","Leading Jet eta",30,-1.5,1.5);
1345 fhLeadingEta->SetXTitle("eta");
1346 fhLeadingEta->SetYTitle("dN/deta");
1347 fhLeadingEta->Sumw2();
1348 fListOfHistos->Add(fhLeadingEta);
1349
1350 fhLeadingPhi = new TH1F("hLeadingPhi","Leading Jet phi",64,-3.2,3.2);
1351 fhLeadingPhi->SetXTitle("phi");
1352 fhLeadingPhi->SetYTitle("dN/dphi");
1353 fhLeadingPhi->Sumw2();
1354 fListOfHistos->Add(fhLeadingPhi);
1355
1356 fhLeadingPt = new TH1F("hLeadingPt","leading Jet p_{T}",50,0,50);
1357 fhLeadingPt->SetXTitle("p_{T} (GeV/c)");
1358 fhLeadingPt->SetYTitle("dN/dp_{T} (1/GeV)");
1359 fhLeadingPt->Sumw2();
1360 fListOfHistos->Add( fhLeadingPt );
1361
1362 fhLeadingPtDiffr = new TH1F("hLeadingPtDiffr","leading Jet p_{T}",50,0,50);
1363 fhLeadingPtDiffr->SetXTitle("P{T} (GeV/c)");
1364 fhLeadingPtDiffr->SetYTitle("dN/dp_{T} (1/GeV)");
1365 fhLeadingPtDiffr->Sumw2();
1366 fListOfHistos->Add( fhLeadingPtDiffr );
1367
1368 fhLeadingEtaMC = new TH1F("hLeadingEtaMC","Leading Jet eta",30,-1.5,1.5);
1369 fhLeadingEtaMC->SetXTitle("eta");
1370 fhLeadingEtaMC->SetYTitle("dN/deta");
1371 fhLeadingEtaMC->Sumw2();
1372 fListOfHistos->Add(fhLeadingEtaMC);
1373
1374 fhLeadingPhiMC = new TH1F("hLeadingPhiMC","Leading Jet phi",64,-3.2,3.2);
1375 fhLeadingPhiMC->SetXTitle("phi");
1376 fhLeadingPhiMC->SetYTitle("dN/dphi");
1377 fhLeadingPhiMC->Sumw2();
1378 fListOfHistos->Add(fhLeadingPhiMC);
1379
1380 fhLeadingPtMC = new TH1F("hLeadingPtMC","Leading Jet p_{T}",50,0,50);
1381 fhLeadingPtMC->SetXTitle("p_{T} (GeV/c)");
1382 fhLeadingPtMC->SetYTitle("dN/dp_{T} (1/GeV)");
1383 fhLeadingPtMC->Sumw2();
1384 fListOfHistos->Add( fhLeadingPtMC );
1385
1386 fhLeadingPtMCDiffr = new TH1F("hLeadingPtMCDiffr","Leading Jet p_{T}",50,0,50);
1387 fhLeadingPtMCDiffr->SetXTitle("p_{T} (GeV/c)");
1388 fhLeadingPtMCDiffr->SetYTitle("dN/dp_{T} (1/GeV)");
1389 fhLeadingPtMCDiffr->Sumw2();
1390 fListOfHistos->Add( fhLeadingPtMCDiffr );
1391
1392 fhEtaTracks = new TH1F("hEtaTracks","eta tracks",30,-1.5,1.5);
1393 fhEtaTracks->SetXTitle("eta");
1394 fhEtaTracks->SetYTitle("dN/deta");
1395 fhEtaTracks->Sumw2();
1396 fListOfHistos->Add(fhEtaTracks);
1397
1398 fhPhiTracks = new TH1F("hPhiTracks","phi tracks",64,-3.2,3.2);
1399 fhPhiTracks->SetXTitle("phi");
1400 fhPhiTracks->SetYTitle("dN/dphi");
1401 fhPhiTracks->Sumw2();
1402 fListOfHistos->Add(fhPhiTracks);
1403
1404 fhPtTracks = new TH1F("hPtTracks","p_{T} tracks",50,0,50);
1405 fhPtTracks->SetXTitle("P{T} (GeV)");
1406 fhPtTracks->SetYTitle("dN/dp_{T} (1/GeV)");
1407 fhPtTracks->Sumw2();
1408 fListOfHistos->Add(fhPtTracks);
1409
1410 fhEtaMCTracks = new TH1F("hEtaMCTracks","eta tracks",30,-1.5,1.5);
1411 fhEtaMCTracks->SetXTitle("eta");
1412 fhEtaMCTracks->SetYTitle("dN/deta");
1413 fhEtaMCTracks->Sumw2();
1414 fListOfHistos->Add(fhEtaMCTracks);
1415
1416 fhPhiMCTracks = new TH1F("hPhiMCTracks","phi tracks",64,-3.2,3.2);
1417 fhPhiMCTracks->SetXTitle("phi");
1418 fhPhiMCTracks->SetYTitle("dN/dphi");
1419 fhPhiMCTracks->Sumw2();
1420 fListOfHistos->Add(fhPhiMCTracks);
1421
1422 fhPtMCTracks = new TH1F("hPtMCTracks","p_{T} tracks",50,0,50);
1423 fhPtMCTracks->SetXTitle("p_{T} (GeV)");
1424 fhPtMCTracks->SetYTitle("dN/dp_{T} (1/GeV)");
1425 fhPtMCTracks->Sumw2();
1426 fListOfHistos->Add(fhPtMCTracks);
1427
1428 fhnTracksVsPtLeading = new TH2F("hnTracksVsPtLeading","",50,0.,50.,20,-0.5,19.5);
1429 fhnTracksVsPtLeading->SetXTitle("p_{T} (GeV/c)");
1430 fhnTracksVsPtLeading->SetYTitle("n tracks");
1431 fhnTracksVsPtLeading->Sumw2();
1432 fListOfHistos->Add( fhnTracksVsPtLeading );
1433
1434 fhdNdEtaPhiDist = new TH1F("hdNdEtaPhiDist","Charge particle density |#eta|< 1 vs #Delta#phi", 120, 0., 2.*TMath::Pi());
1435 fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
1436 fhdNdEtaPhiDist->SetYTitle("dN{ch}/d#etad#phi");
1437 fhdNdEtaPhiDist->Sumw2();
1438 fListOfHistos->Add( fhdNdEtaPhiDist );
1439
1440 fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt","P{T}^{90, max} vs Leading Jet P{T}",50,0.,50.);
1441 fhRegionSumPtMaxVsEt->SetXTitle("P{T} (GeV/c)");
1442 fhRegionSumPtMaxVsEt->Sumw2();
1443 fListOfHistos->Add( fhRegionSumPtMaxVsEt );
1444
1445 fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt","N{ch}^{90, max} vs Leading Jet P{T}",50,0.,50.);
1446 fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
1447 fhRegionMultMaxVsEt->Sumw2();
1448 fListOfHistos->Add( fhRegionMultMaxVsEt );
1449
1450 fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt","P{T}^{90, min} vs Leading Jet P{T}",50,0.,50.);
1451 fhRegionSumPtMinVsEt->SetXTitle("P{T} (GeV/c)");
1452 fhRegionSumPtMinVsEt->Sumw2();
1453 fListOfHistos->Add( fhRegionSumPtMinVsEt );
1454
1455 fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt","N{ch}^{90, min} vs Leading Jet P{T}",50,0.,50.);
1456 fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
1457 fhRegionMultMinVsEt->Sumw2();
1458 fListOfHistos->Add( fhRegionMultMinVsEt );
1459
1460 // V0s
1461
1462 fhNV0s = new TH1F("hNV0s","n V0s",50,0,50);
1463 fhNV0s->SetXTitle("n V0s");
1464 fhNV0s->Sumw2();
1465 fListOfHistos->Add(fhNV0s);
1466
1467 fhV0onFly = new TH1F("hV0onFly","on-the-fly V0",5,0,5);
1468 fhV0onFly->SetXTitle("is on-the-fly V0");
1469 fhV0onFly->Sumw2();
1470 fListOfHistos->Add(fhV0onFly);
1471
1472 fhV0DCADaughters = new TH1F("hV0DCADaughters","V0 DCA daughters",200,0,2.0);
1473 fhV0DCADaughters->SetXTitle("V0 DCA daughters");
1474 fhV0DCADaughters->Sumw2();
1475 fListOfHistos->Add(fhV0DCADaughters);
1476
1477 fhV0Radius = new TH1F("hV0Radius","V0 radius",2500,0,250);
1478 fhV0Radius->SetXTitle("V0 radius");
1479 fhV0Radius->Sumw2();
1480 fListOfHistos->Add(fhV0Radius);
1481
1482 fhV0DCAToVertex = new TH1F("hV0DCAToVertex","",100,0,10);
1483 fhV0DCAToVertex->SetXTitle("V0 DCA (cm)");
1484 fhV0DCAToVertex->Sumw2();
1485 fListOfHistos->Add(fhV0DCAToVertex);
1486
1487 fhV0DCAToVertexK0 = new TH1F("hV0DCAToVertexK0","",100,0,10);
1488 fhV0DCAToVertexK0->SetXTitle("V0 DCA (cm)");
1489 fhV0DCAToVertexK0->Sumw2();
1490 fListOfHistos->Add(fhV0DCAToVertexK0);
1491
1492 fhV0InvMassK0 = new TH1F("hV0InvMassK0","",2000,0,2);
1493 fhV0InvMassK0->SetXTitle("inv mass (GeV)");
1494 fhV0InvMassK0->Sumw2();
1495 fListOfHistos->Add(fhV0InvMassK0);
1496
1497 fhV0InvMassK0JetEvt = new TH1F("hV0InvMassK0JetEvt","",2000,0,2);
1498 fhV0InvMassK0JetEvt->SetXTitle("inv mass (GeV)");
1499 fhV0InvMassK0JetEvt->Sumw2();
1500 fListOfHistos->Add(fhV0InvMassK0JetEvt);
1501
1502 fhV0InvMassLambda = new TH1F("hV0InvMassLambda","",2000,0,2);
1503 fhV0InvMassLambda->SetXTitle("inv mass (GeV)");
1504 fhV0InvMassLambda->Sumw2();
1505 fListOfHistos->Add(fhV0InvMassLambda);
1506
1507 fhV0InvMassAntiLambda = new TH1F("hV0InvMassAntiLambda","",2000,0,2);
1508 fhV0InvMassAntiLambda->SetXTitle("inv mass (GeV)");
1509 fhV0InvMassAntiLambda->Sumw2();
1510 fListOfHistos->Add(fhV0InvMassAntiLambda);
1511
1512 fhV0InvMassLambdaJetEvt = new TH1F("hV0InvMassLambdaJetEvt","",2000,0,2);
1513 fhV0InvMassLambdaJetEvt->SetXTitle("inv mass (GeV)");
1514 fhV0InvMassLambdaJetEvt->Sumw2();
1515 fListOfHistos->Add(fhV0InvMassLambdaJetEvt);
1516
1517 fhV0InvMassAntiLambdaJetEvt = new TH1F("hV0InvMassAntiLambdaJetEvt","",2000,0,2);
1518 fhV0InvMassAntiLambdaJetEvt->SetXTitle("inv mass (GeV)");
1519 fhV0InvMassAntiLambdaJetEvt->Sumw2();
1520 fListOfHistos->Add(fhV0InvMassAntiLambdaJetEvt);
1521
1522 fhdROpanK0VsPt = new TH2F("hdROpanK0VsPt","V0 dR vs pt",100,0,10,100,0,10);
1523 fhdROpanK0VsPt->SetXTitle("opening angle R (rad)");
1524 fhdROpanK0VsPt->Sumw2();
1525 fListOfHistos->Add(fhdROpanK0VsPt);
1526
1527 fhdPhiJetV0 = new TH1F("hdPhiJetV0","",640,-3.2,3.2);
1528 fhdPhiJetV0->SetXTitle("#Delta #phi V0 - jet");
1529 fhdPhiJetV0->Sumw2();
1530 fListOfHistos->Add(fhdPhiJetV0);
1531
1532 fhdPhiJetK0 = new TH1F("hdPhiJetK0","",640,-3.2,3.2);
1533 fhdPhiJetK0->SetXTitle("#Delta #phi V0 - jet");
1534 fhdPhiJetK0->Sumw2();
1535 fListOfHistos->Add(fhdPhiJetK0);
1536
1537 fhdRJetK0 = new TH1F("hdRJetK0","dN/dR K0-jet",500,0,5);
1538 fhdRJetK0->SetXTitle("#Delta R K0 - jet");
1539 fhdRJetK0->SetYTitle("1/N{jet} dN/dR");
1540 fhdRJetK0->Sumw2();
1541 fListOfHistos->Add(fhdRJetK0);
1542
1543 fhdNdptV0 = new TH1F("hdNdptV0","dN/dpt V0",100,0,10);
1544 fhdNdptV0->SetXTitle("p_{T} (GeV/c)");
1545 fhdNdptV0->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1546 fhdNdptV0->Sumw2();
1547 fListOfHistos->Add(fhdNdptV0);
1548
1549 fhdNdptK0 = new TH1F("hdNdptK0","",100,0,10);
1550 fhdNdptK0->SetXTitle("p_{T} (GeV/c)");
1551 fhdNdptK0->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1552 fhdNdptK0->Sumw2();
1553 fListOfHistos->Add(fhdNdptK0);
1554
1555 fhPtVsEtaK0 = new TH2F("hPtVsEtaK0","",20,-1,1,100,0,10);
1556 fhPtVsEtaK0->SetXTitle("#eta");
1557 fhPtVsEtaK0->SetYTitle("p_{T} (GeV/c)");
1558 fhPtVsEtaK0->Sumw2();
1559 fListOfHistos->Add(fhPtVsEtaK0);
1560
1561 fhV0InvMassK0DCA = new TH1F("hV0InvMassK0DCA","",2000,0,2);
1562 fhV0InvMassK0DCA->SetXTitle("inv mass (GeV)");
1563 fhV0InvMassK0DCA->Sumw2();
1564 fListOfHistos->Add(fhV0InvMassK0DCA);
1565
1566 fhV0InvMassK0DCAdEdx = new TH1F("hV0InvMassK0DCAdEdx","",2000,0,2);
1567 fhV0InvMassK0DCAdEdx->SetXTitle("inv mass (GeV)");
1568 fhV0InvMassK0DCAdEdx->Sumw2();
1569 fListOfHistos->Add(fhV0InvMassK0DCAdEdx);
1570
1571 fhV0InvMassK0DCAPID = new TH1F("hV0InvMassK0DCAPID","",2000,0,2);
1572 fhV0InvMassK0DCAPID->SetXTitle("inv mass (GeV)");
1573 fhV0InvMassK0DCAPID->Sumw2();
1574 fListOfHistos->Add(fhV0InvMassK0DCAPID);
1575
1576 fhV0InvMassLambdaDCAdEdx = new TH1F("hV0InvMassLambdaDCAdEdx","",2000,0,2);
1577 fhV0InvMassLambdaDCAdEdx->SetXTitle("inv mass (GeV)");
1578 fhV0InvMassLambdaDCAdEdx->Sumw2();
1579 fListOfHistos->Add(fhV0InvMassLambdaDCAdEdx);
1580
1581 fhV0InvMassAntiLambdaDCAdEdx = new TH1F("hV0InvMassAntiLambdaDCAdEdx","",2000,0,2);
1582 fhV0InvMassAntiLambdaDCAdEdx->SetXTitle("inv mass (GeV)");
1583 fhV0InvMassAntiLambdaDCAdEdx->Sumw2();
1584 fListOfHistos->Add(fhV0InvMassAntiLambdaDCAdEdx);
1585
1586 fhdNdptK0DCA = new TH1F("hdNdptK0DCA","",100,0,10);
1587 fhdNdptK0DCA->SetXTitle("p_{T} (GeV/c)");
1588 fhdNdptK0DCA->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1589 fhdNdptK0DCA->Sumw2();
1590 fListOfHistos->Add(fhdNdptK0DCA);
1591
1592 fhdNdptK0DCAdEdx = new TH1F("hdNdptK0DCAdEdx","",100,0,10);
1593 fhdNdptK0DCAdEdx->SetXTitle("p_{T} (GeV/c)");
1594 fhdNdptK0DCAdEdx->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1595 fhdNdptK0DCAdEdx->Sumw2();
1596 fListOfHistos->Add(fhdNdptK0DCAdEdx);
1597
1598 fhV0InvMassK0Min = new TH1F("hV0InvMassK0Min","",2000,0,2);
1599 fhV0InvMassK0Min->SetXTitle("inv mass (GeV)");
1600 fhV0InvMassK0Min->Sumw2();
1601 fListOfHistos->Add(fhV0InvMassK0Min);
1602
1603 fhV0InvMassLambdaMin = new TH1F("hV0InvMassLambdaMin","",2000,0,2);
1604 fhV0InvMassLambdaMin->SetXTitle("inv mass (GeV)");
1605 fhV0InvMassLambdaMin->Sumw2();
1606 fListOfHistos->Add(fhV0InvMassLambdaMin);
1607
1608 fhV0InvMassAntiLambdaMin = new TH1F("hV0InvMassAntiLambdaMin","",2000,0,2);
1609 fhV0InvMassAntiLambdaMin->SetXTitle("inv mass (GeV)");
1610 fhV0InvMassAntiLambdaMin->Sumw2();
1611 fListOfHistos->Add(fhV0InvMassAntiLambdaMin);
1612
1613 fhV0InvMassK0Max = new TH1F("hV0InvMassK0Max","",2000,0,2);
1614 fhV0InvMassK0Max->SetXTitle("inv mass (GeV)");
1615 fhV0InvMassK0Max->Sumw2();
1616 fListOfHistos->Add(fhV0InvMassK0Max);
1617
1618 fhV0InvMassLambdaMax = new TH1F("hV0InvMassLambdaMax","",2000,0,2);
1619 fhV0InvMassLambdaMax->SetXTitle("inv mass (GeV)");
1620 fhV0InvMassLambdaMax->Sumw2();
1621 fListOfHistos->Add(fhV0InvMassLambdaMax);
1622
1623 fhV0InvMassAntiLambdaMax = new TH1F("hV0InvMassAntiLambdaMax","",2000,0,2);
1624 fhV0InvMassAntiLambdaMax->SetXTitle("inv mass (GeV)");
1625 fhV0InvMassAntiLambdaMax->Sumw2();
1626 fListOfHistos->Add(fhV0InvMassAntiLambdaMax);
1627
1628 fhV0InvMassK0Jet = new TH1F("hV0InvMassK0Jet","",2000,0,2);
1629 fhV0InvMassK0Jet->SetXTitle("inv mass (GeV)");
1630 fhV0InvMassK0Jet->Sumw2();
1631 fListOfHistos->Add(fhV0InvMassK0Jet);
1632
1633 fhV0InvMassLambdaJet = new TH1F("hV0InvMassLambdaJet","",2000,0,2);
1634 fhV0InvMassLambdaJet->SetXTitle("inv mass (GeV)");
1635 fhV0InvMassLambdaJet->Sumw2();
1636 fListOfHistos->Add(fhV0InvMassLambdaJet);
1637
1638 fhV0InvMassAntiLambdaJet = new TH1F("hV0InvMassAntiLambdaJet","",2000,0,2);
1639 fhV0InvMassAntiLambdaJet->SetXTitle("inv mass (GeV)");
1640 fhV0InvMassAntiLambdaJet->Sumw2();
1641 fListOfHistos->Add(fhV0InvMassAntiLambdaJet);
1642
1643 fhV0InvMassK0Lambda = new TH1F("hV0InvMassK0Lambda","",2000,0,2);
1644 fhV0InvMassK0Lambda->SetXTitle("inv mass (GeV)");
1645 fhV0InvMassK0Lambda->Sumw2();
1646 fListOfHistos->Add(fhV0InvMassK0Lambda);
1647
1648 fhdNdptK0JetEvt = new TH1F("hdNdptK0JetEvt","",100,0,10);
1649 fhdNdptK0JetEvt->SetXTitle("p_{T} (GeV/c)");
1650 fhdNdptK0JetEvt->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1651 fhdNdptK0JetEvt->Sumw2();
1652 fListOfHistos->Add(fhdNdptK0JetEvt);
1653
1654 fhdNdzK0 = new TH1F("hdNdzK0","",150,0,1.5);
1655 fhdNdzK0->SetXTitle("z");
1656 fhdNdzK0->SetYTitle("1/N{jet} dN/dz");
1657 fhdNdzK0->Sumw2();
1658 fListOfHistos->Add(fhdNdzK0);
1659
1660 fhdNdxiK0 = new TH1F("hdNdxiK0","",100,0,10);
1661 fhdNdxiK0->SetXTitle("xi");
1662 fhdNdxiK0->SetYTitle("1/N{jet} dN/dxi");
1663 fhdNdxiK0->Sumw2();
1664 fListOfHistos->Add(fhdNdxiK0);
1665
1666 fhdNdzK0Max = new TH1F("hdNdzK0Max","",150,0,1.5);
1667 fhdNdzK0Max->SetXTitle("z");
1668 fhdNdzK0Max->SetYTitle("1/N{jet} dN/dz");
1669 fhdNdzK0Max->Sumw2();
1670 fListOfHistos->Add(fhdNdzK0Max);
1671
1672 fhdNdxiK0Max = new TH1F("hdNdxiK0Max","",100,0,10);
1673 fhdNdxiK0Max->SetXTitle("xi");
1674 fhdNdxiK0Max->SetYTitle("1/N{jet} dN/dxi");
1675 fhdNdxiK0Max->Sumw2();
1676 fListOfHistos->Add(fhdNdxiK0Max);
1677
1678 fhdNdzLambdaMax = new TH1F("hdNdzLambdaMax","",150,0,1.5);
1679 fhdNdzLambdaMax->SetXTitle("z");
1680 fhdNdzLambdaMax->SetYTitle("1/N{jet} dN/dz");
1681 fhdNdzLambdaMax->Sumw2();
1682 fListOfHistos->Add(fhdNdzLambdaMax);
1683
1684 fhdNdxiLambdaMax = new TH1F("hdNdxiLambdaMax","",700,0,7);
1685 fhdNdxiLambdaMax->SetXTitle("xi");
1686 fhdNdxiLambdaMax->SetYTitle("1/N{jet} dN/dxi");
1687 fhdNdxiLambdaMax->Sumw2();
1688 fListOfHistos->Add(fhdNdxiLambdaMax);
1689
1690 fhdNdptK0Max = new TH1F("hdNdptK0Max","",100,0,10);
1691 fhdNdptK0Max->SetXTitle("p_{T} (GeV/c)");
1692 fhdNdptK0Max->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1693 fhdNdptK0Max->Sumw2();
1694 fListOfHistos->Add(fhdNdptK0Max);
1695
1696 fhdNdptLambdaMax = new TH1F("hdNdptLambdaMax","",100,0,10);
1697 fhdNdptLambdaMax->SetXTitle("p_{T} (GeV/c)");
1698 fhdNdptLambdaMax->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1699 fhdNdptLambdaMax->Sumw2();
1700 fListOfHistos->Add(fhdNdptLambdaMax);
1701
1702 fhdNdzK0Min = new TH1F("hdNdzK0Min","",150,0,1.5);
1703 fhdNdzK0Min->SetXTitle("z");
1704 fhdNdzK0Min->SetYTitle("1/N{jet} dN/dz");
1705 fhdNdzK0Min->Sumw2();
1706 fListOfHistos->Add(fhdNdzK0Min);
1707
1708 fhdNdxiK0Min = new TH1F("hdNdxiK0Min","",100,0,10);
1709 fhdNdxiK0Min->SetXTitle("xi");
1710 fhdNdxiK0Min->SetYTitle("1/N{jet} dN/dxi");
1711 fhdNdxiK0Min->Sumw2();
1712 fListOfHistos->Add(fhdNdxiK0Min);
1713
1714 fhdNdzLambdaMin = new TH1F("hdNdzLambdaMin","",150,0,1.5);
1715 fhdNdzLambdaMin->SetXTitle("z");
1716 fhdNdzLambdaMin->SetYTitle("1/N{jet} dN/dz");
1717 fhdNdzLambdaMin->Sumw2();
1718 fListOfHistos->Add(fhdNdzLambdaMin);
1719
1720 fhdNdxiLambdaMin = new TH1F("hdNdxiLambdaMin","",700,0,7);
1721 fhdNdxiLambdaMin->SetXTitle("xi");
1722 fhdNdxiLambdaMin->SetYTitle("1/N{jet} dN/dxi");
1723 fhdNdxiLambdaMin->Sumw2();
1724 fListOfHistos->Add(fhdNdxiLambdaMin);
1725
1726 fhdNdptK0Min = new TH1F("hdNdptK0Min","",100,0,10);
1727 fhdNdptK0Min->SetXTitle("p_{T} (GeV/c)");
1728 fhdNdptK0Min->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1729 fhdNdptK0Min->Sumw2();
1730 fListOfHistos->Add(fhdNdptK0Min);
1731
1732 fhdNdptLambdaMin = new TH1F("hdNdptLambdaMin","",100,0,10);
1733 fhdNdptLambdaMin->SetXTitle("p_{T} (GeV/c)");
1734 fhdNdptLambdaMin->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1735 fhdNdptLambdaMin->Sumw2();
1736 fListOfHistos->Add(fhdNdptLambdaMin);
1737
1738 fhdNdzK0Jet = new TH1F("hdNdzK0Jet","",150,0,1.5);
1739 fhdNdzK0Jet->SetXTitle("z");
1740 fhdNdzK0Jet->SetYTitle("1/N{jet} dN/dz");
1741 fhdNdzK0Jet->Sumw2();
1742 fListOfHistos->Add(fhdNdzK0Jet);
1743
1744 fhdNdxiK0Jet = new TH1F("hdNdxiK0Jet","",100,0,10);
1745 fhdNdxiK0Jet->SetXTitle("xi");
1746 fhdNdxiK0Jet->SetYTitle("1/N{jet} dN/dxi");
1747 fhdNdxiK0Jet->Sumw2();
1748 fListOfHistos->Add(fhdNdxiK0Jet);
1749
1750 fhdNdzLambdaJet = new TH1F("hdNdzLambdaJet","",150,0,1.5);
1751 fhdNdzLambdaJet->SetXTitle("z");
1752 fhdNdzLambdaJet->SetYTitle("1/N{jet} dN/dz");
1753 fhdNdzLambdaJet->Sumw2();
1754 fListOfHistos->Add(fhdNdzLambdaJet);
1755
1756 fhdNdxiLambdaJet = new TH1F("hdNdxiLambdaJet","",700,0,7);
1757 fhdNdxiLambdaJet->SetXTitle("xi");
1758 fhdNdxiLambdaJet->SetYTitle("1/N{jet} dN/dxi");
1759 fhdNdxiLambdaJet->Sumw2();
1760 fListOfHistos->Add(fhdNdxiLambdaJet);
1761
1762 fhdNdptK0Jet = new TH1F("hdNdptK0Jet","",100,0,10);
1763 fhdNdptK0Jet->SetXTitle("p_{T} (GeV/c)");
1764 fhdNdptK0Jet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1765 fhdNdptK0Jet->Sumw2();
1766 fListOfHistos->Add(fhdNdptK0Jet);
1767
1768 fhdNdptLambdaJet = new TH1F("hdNdptLambdaJet","",100,0,10);
1769 fhdNdptLambdaJet->SetXTitle("p_{T} (GeV/c)");
1770 fhdNdptLambdaJet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1771 fhdNdptLambdaJet->Sumw2();
1772 fListOfHistos->Add(fhdNdptLambdaJet);
1773
1774 fhdEdxVsMomV0 = new TH2F("hdEdxVsMomV0","",200,-10,10,200,0,2000);
1775 fhdEdxVsMomV0->SetXTitle("mom (GeV/c)");
1776 fhdEdxVsMomV0->SetYTitle("dE/dx");
1777 fhdEdxVsMomV0->Sumw2();
1778 fListOfHistos->Add(fhdEdxVsMomV0);
1779
1780 fhdEdxVsMomV0pidEdx = new TH2F("hdEdxVsMomV0pidEdx","",200,-10,10,200,0,2000);
1781 fhdEdxVsMomV0pidEdx->SetXTitle("mom (GeV/c)");
1782 fhdEdxVsMomV0pidEdx->SetYTitle("dE/dx");
1783 fhdEdxVsMomV0pidEdx->Sumw2();
1784 fListOfHistos->Add(fhdEdxVsMomV0pidEdx);
1785
1786 fhdEdxVsMomV0piPID = new TH2F("hdEdxVsMomV0piPID","",200,-10,10,200,0,2000);
1787 fhdEdxVsMomV0piPID->SetXTitle("mom (GeV/c)");
1788 fhdEdxVsMomV0piPID->SetYTitle("dE/dx");
1789 fhdEdxVsMomV0piPID->Sumw2();
1790 fListOfHistos->Add(fhdEdxVsMomV0piPID);
1791
1792 fhdPhiJetK0MC = new TH1F("hdPhiJetK0MC","",640,-3.2,3.2);
1793 fhdPhiJetK0MC->SetXTitle("#Delta #phi K0 - jet");
1794 fhdPhiJetK0MC->SetYTitle("1/N{jet} dN/dphi");
1795 fhdPhiJetK0MC->Sumw2();
1796 fListOfHistos->Add(fhdPhiJetK0MC);
1797
1798 fhdRJetK0MC = new TH1F("hdRJetK0MC","dN/R K0-jet",500,0,5);
1799 fhdRJetK0MC->SetXTitle("#Delta R K0 - jet");
1800 fhdRJetK0MC->SetYTitle("1/N{jet} dN/dR");
1801 fhdRJetK0MC->Sumw2();
1802 fListOfHistos->Add(fhdRJetK0MC);
1803
1804 fhdNdptK0MCMax = new TH1F("hdNdptK0MCMax","",100,0,10);
1805 fhdNdptK0MCMax->SetXTitle("p_{T} (GeV/c)");
1806 fhdNdptK0MCMax->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1807 fhdNdptK0MCMax->Sumw2();
1808 fListOfHistos->Add(fhdNdptK0MCMax);
1809
1810 fhdNdptLambdaLambdaBarMCMax = new TH1F("hdNdptLambdaLambdaBarMCMax","",100,0,10);
1811 fhdNdptLambdaLambdaBarMCMax->SetXTitle("p_{T} (GeV/c)");
1812 fhdNdptLambdaLambdaBarMCMax->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1813 fhdNdptLambdaLambdaBarMCMax->Sumw2();
1814 fListOfHistos->Add(fhdNdptLambdaLambdaBarMCMax);
1815
1816 fhdNdptK0MCMin = new TH1F("hdNdptK0MCMin","",100,0,10);
1817 fhdNdptK0MCMin->SetXTitle("p_{T} (GeV/c)");
1818 fhdNdptK0MCMin->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1819 fhdNdptK0MCMin->Sumw2();
1820 fListOfHistos->Add(fhdNdptK0MCMin);
1821
1822 fhdNdptLambdaLambdaBarMCMin = new TH1F("hdNdptLambdaLambdaBarMCMin","",100,0,10);
1823 fhdNdptLambdaLambdaBarMCMin->SetXTitle("p_{T} (GeV/c)");
1824 fhdNdptLambdaLambdaBarMCMin->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1825 fhdNdptLambdaLambdaBarMCMin->Sumw2();
1826 fListOfHistos->Add(fhdNdptLambdaLambdaBarMCMin);
1827
1828 fhdNdptOmegaMCMin = new TH1F("hdNdptOmegaMCMin","",100,0,10);;
1829 fhdNdptOmegaMCMin->SetXTitle("p_{T} (GeV/c)");
1830 fhdNdptOmegaMCMin->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1831 fhdNdptOmegaMCMin->Sumw2();
1832 fListOfHistos->Add(fhdNdptOmegaMCMin);
1833
1834 fhdNdptOmegaBarMCMin = new TH1F("hdNdptOmegaBarMCMin","",100,0,10);;
1835 fhdNdptOmegaBarMCMin->SetXTitle("p_{T} (GeV/c)");
1836 fhdNdptOmegaBarMCMin->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1837 fhdNdptOmegaBarMCMin->Sumw2();
1838 fListOfHistos->Add(fhdNdptOmegaBarMCMin);
1839
1840 fhdNdptchPiMCMin = new TH1F("hdNdptchPiMCMin","",100,0,10);
1841 fhdNdptchPiMCMin->SetXTitle("p_{T} (GeV/c)");
1842 fhdNdptchPiMCMin->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1843 fhdNdptchPiMCMin->Sumw2();
1844 fListOfHistos->Add(fhdNdptchPiMCMin);
1845
1846 fhdNdptchKMCMin = new TH1F("hdNdptchKMCMin","",100,0,10);
1847 fhdNdptchKMCMin->SetXTitle("p_{T} (GeV/c)");
1848 fhdNdptchKMCMin->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1849 fhdNdptchKMCMin->Sumw2();
1850 fListOfHistos->Add(fhdNdptchKMCMin);
1851
1852 fhdNdptpMCMin = new TH1F("hdNdptpMCMin","",100,0,10);
1853 fhdNdptpMCMin->SetXTitle("p_{T} (GeV/c)");
1854 fhdNdptpMCMin->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1855 fhdNdptpMCMin->Sumw2();
1856 fListOfHistos->Add(fhdNdptpMCMin);
1857
1858 fhdNdptpBarMCMin = new TH1F("hdNdptpBarMCMin","",100,0,10);
1859 fhdNdptpBarMCMin->SetXTitle("p_{T} (GeV/c)");
1860 fhdNdptpBarMCMin->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1861 fhdNdptpBarMCMin->Sumw2();
1862 fListOfHistos->Add(fhdNdptpBarMCMin);
1863
1864 fhdNdptK0MCJet = new TH1F("hdNdptK0MCJet","",100,0,10);
1865 fhdNdptK0MCJet->SetXTitle("p_{T} (GeV/c)");
1866 fhdNdptK0MCJet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1867 fhdNdptK0MCJet->Sumw2();
1868 fListOfHistos->Add(fhdNdptK0MCJet);
1869
1870 fhdNdptLambdaLambdaBarMCJet = new TH1F("hdNdptLambdaLambdaBarMCJet","",100,0,10);
1871 fhdNdptLambdaLambdaBarMCJet->SetXTitle("p_{T} (GeV/c)");
1872 fhdNdptLambdaLambdaBarMCJet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1873 fhdNdptLambdaLambdaBarMCJet->Sumw2();
1874 fListOfHistos->Add(fhdNdptLambdaLambdaBarMCJet);
1875
1876 fhdNdptchPiMCJet = new TH1F("hdNdptchPiMCJet","",100,0,10);
1877 fhdNdptchPiMCJet->SetXTitle("p_{T} (GeV/c)");
1878 fhdNdptchPiMCJet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1879 fhdNdptchPiMCJet->Sumw2();
1880 fListOfHistos->Add(fhdNdptchPiMCJet);
1881
1882 fhdNdptchKMCJet = new TH1F("hdNdptchKMCJet","",100,0,10);
1883 fhdNdptchKMCJet->SetXTitle("p_{T} (GeV/c)");
1884 fhdNdptchKMCJet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1885 fhdNdptchKMCJet->Sumw2();
1886 fListOfHistos->Add(fhdNdptchKMCJet);
1887
1888 fhdNdptpMCJet = new TH1F("hdNdptpMCJet","",100,0,10);
1889 fhdNdptpMCJet->SetXTitle("p_{T} (GeV/c)");
1890 fhdNdptpMCJet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1891 fhdNdptpMCJet->Sumw2();
1892 fListOfHistos->Add(fhdNdptpMCJet);
1893
1894 fhdNdptpBarMCJet = new TH1F("hdNdptpBarMCJet","",100,0,10);
1895 fhdNdptpBarMCJet->SetXTitle("p_{T} (GeV/c)");
1896 fhdNdptpBarMCJet->SetYTitle("1/N{jet} dN/dpt (1/GeV)");
1897 fhdNdptpBarMCJet->Sumw2();
1898 fListOfHistos->Add(fhdNdptpBarMCJet);
1899
1900 fhdRV0MC = new TH1F("hdRV0MC","",500,0.,1.);
1901 fhdRV0MC->SetXTitle("#Delta R");
1902 fhdRV0MC->SetYTitle("");
1903 fhdRV0MC->Sumw2();
1904 fListOfHistos->Add(fhdRV0MC);
1905
1906
1907 fhPIDMC = CreatePIDhisto("hPIDMC");
1908 fhPIDMC->Sumw2();
1909 fListOfHistos->Add(fhPIDMC);
1910
1911 fhPIDMCAll = CreatePIDhisto("hPIDMCAll");
1912 fhPIDMCAll->Sumw2();
1913 fListOfHistos->Add(fhPIDMCAll);
1914
1915 fhPIDMCMin = CreatePIDhisto("hPIDMCMin");
1916 fhPIDMCMin->Sumw2();
1917 fListOfHistos->Add(fhPIDMCMin);
1918
1919 fhPIDMCJet = CreatePIDhisto("hPIDMCJet");
1920 fhPIDMCJet->Sumw2();
1921 fListOfHistos->Add(fhPIDMCJet);
1922
1923 fhdNdptgammaMC = new TH1F("hdNdptgammaMC","",100,0,10);
1924 fhdNdptgammaMC->SetXTitle("p_{T} (GeV/c)");
1925 fhdNdptgammaMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1926 fhdNdptgammaMC->Sumw2();
1927 fListOfHistos->Add(fhdNdptgammaMC);
1928
1929 fhdNdptchPiMC = new TH1F("hdNdptchPiMC","",100,0,10);;
1930 fhdNdptchPiMC->SetXTitle("p_{T} (GeV/c)");
1931 fhdNdptchPiMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1932 fhdNdptchPiMC->Sumw2();
1933 fListOfHistos->Add(fhdNdptchPiMC);
1934
1935 fhdNdptpi0MC = new TH1F("hdNdptpi0MC","",100,0,10);;
1936 fhdNdptpi0MC->SetXTitle("p_{T} (GeV/c)");
1937 fhdNdptpi0MC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1938 fhdNdptpi0MC->Sumw2();
1939 fListOfHistos->Add(fhdNdptpi0MC);
1940
1941 fhdNdptK0MC = new TH1F("hdNdptK0MC","",100,0,10);;
1942 fhdNdptK0MC->SetXTitle("p_{T} (GeV/c)");
1943 fhdNdptK0MC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1944 fhdNdptK0MC->Sumw2();
1945 fListOfHistos->Add(fhdNdptK0MC);
1946
1947 fhdNdptchKMC = new TH1F("hdNdptchKMC","",100,0,10);;
1948 fhdNdptchKMC->SetXTitle("p_{T} (GeV/c)");
1949 fhdNdptchKMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1950 fhdNdptchKMC->Sumw2();
1951 fListOfHistos->Add(fhdNdptchKMC);
1952
1953 fhdNdptpMC = new TH1F("hdNdptpMC","",100,0,10);;
1954 fhdNdptpMC->SetXTitle("p_{T} (GeV/c)");
1955 fhdNdptpMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1956 fhdNdptpMC->Sumw2();
1957 fListOfHistos->Add(fhdNdptpMC);
1958
1959 fhdNdptpBarMC = new TH1F("hdNdptpBarMC","",100,0,10);;
1960 fhdNdptpBarMC->SetXTitle("p_{T} (GeV/c)");
1961 fhdNdptpBarMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1962 fhdNdptpBarMC->Sumw2();
1963 fListOfHistos->Add(fhdNdptpBarMC);
1964
1965 fhdNdptLambdaLambdaBarMC = new TH1F("hdNdptLambdaLambdaBarMC","",100,0,10);;
1966 fhdNdptLambdaLambdaBarMC->SetXTitle("p_{T} (GeV/c)");
1967 fhdNdptLambdaLambdaBarMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1968 fhdNdptLambdaLambdaBarMC->Sumw2();
1969 fListOfHistos->Add(fhdNdptLambdaLambdaBarMC);
1970
1971 fhdNdptOmegaMC = new TH1F("hdNdptOmegaMC","",100,0,10);;
1972 fhdNdptOmegaMC->SetXTitle("p_{T} (GeV/c)");
1973 fhdNdptOmegaMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1974 fhdNdptOmegaMC->Sumw2();
1975 fListOfHistos->Add(fhdNdptOmegaMC);
1976
1977 fhdNdptOmegaBarMC = new TH1F("hdNdptOmegaBarMC","",100,0,10);;
1978 fhdNdptOmegaBarMC->SetXTitle("p_{T} (GeV/c)");
1979 fhdNdptOmegaBarMC->SetYTitle("1/N{event} dN/dpt (1/GeV)");
1980 fhdNdptOmegaBarMC->Sumw2();
1981 fListOfHistos->Add(fhdNdptOmegaBarMC);
1982
1983 fhdNdxiMC = new TH1F("hdNdxiMC","",100,0,10);
1984 fhdNdxiMC->SetXTitle("#xi");
1985 fhdNdxiMC->Sumw2();
1986 fListOfHistos->Add(fhdNdxiMC);
1987
1988 fhdNdxiK0MC = new TH1F("hdNdxiK0MC","",100,0,10);
1989 fhdNdxiK0MC->SetXTitle("#xi");
1990 fhdNdxiK0MC->Sumw2();
1991 fListOfHistos->Add(fhdNdxiK0MC);
1992
1993 fhdNdxiK0MCJet = new TH1F("hdNdxiK0MCJet","",100,0,10);
1994 fhdNdxiK0MCJet->SetXTitle("#xi");
1995 fhdNdxiK0MCJet->Sumw2();
1996 fListOfHistos->Add(fhdNdxiK0MCJet);
1997
1998 fhdNdzK0MC = new TH1F("hdNdzK0MC","",150,0,1.5);
1999 fhdNdzK0MC->SetXTitle("z");
2000 fhdNdzK0MC->SetYTitle("1/N{jet} dN/dz");
2001 fhdNdzK0MC->Sumw2();
2002 fListOfHistos->Add(fhdNdzK0MC);
2003
2004 fhdNdzK0MCJet = new TH1F("hdNdzK0MCJet","",150,0,1.5);
2005 fhdNdzK0MCJet->SetXTitle("z");
2006 fhdNdzK0MCJet->SetYTitle("1/N{jet} dN/dz");
2007 fhdNdzK0MCJet->Sumw2();
2008 fListOfHistos->Add(fhdNdzK0MCJet);
2009
2010 fhdNdptK0MCJetEvt = new TH1F("hdNdptK0MCJetEvt","",100,0,10);;
2011 fhdNdptK0MCJetEvt->SetXTitle("p_{T} (GeV/c)");
2012 fhdNdptK0MCJetEvt->SetYTitle("1/N{event} dN/dpt (1/GeV)");
2013 fhdNdptK0MCJetEvt->Sumw2();
2014 fListOfHistos->Add(fhdNdptK0MCJetEvt);
2015
2016 fhnJetsAODvsMC = new TH2F("hnJetsAODvsMC","",20,0,20,20,0.,20.);
2017 fhnJetsAODvsMC->SetXTitle("n jets MC");
2018 fhnJetsAODvsMC->SetXTitle("n jets AOD");
2019 fhnJetsAODvsMC->Sumw2();
2020 fListOfHistos->Add(fhnJetsAODvsMC);
2021
2022 fhLeadingPtAODvsMC = new TH2F("hLeadingPtAODvsMC","",20,0,20,20,0.,20.);
2023 fhLeadingPtAODvsMC->SetXTitle("p_{T} MC (GeV/c)");
2024 fhLeadingPtAODvsMC->SetYTitle("p_{T} AOD (GeV/c)");
2025 fhLeadingPtAODvsMC->Sumw2();
2026 fListOfHistos->Add(fhLeadingPtAODvsMC);
2027
2028 fhLeadingEtaAODvsMC = new TH2F("hLeadingEtaAODvsMC","",20,-1,1,20,-1.,1.);
2029 fhLeadingEtaAODvsMC->SetXTitle("#eta MC");
2030 fhLeadingEtaAODvsMC->SetYTitle("#eta AOD");
2031 fhLeadingEtaAODvsMC->Sumw2();
2032 fListOfHistos->Add(fhLeadingEtaAODvsMC);
2033
2034 fhLeadingPhiAODvsMC = new TH2F("hLeadingPhiAODvsMC","",63,0,6.3,63,0.,6.3);
2035 fhLeadingPhiAODvsMC->SetXTitle("#phi MC");
2036 fhLeadingPhiAODvsMC->SetYTitle("#phi AOD");
2037 fhLeadingPhiAODvsMC->Sumw2();
2038 fListOfHistos->Add(fhLeadingPhiAODvsMC);
2039
2040
2041 fhnTracksLeadingAODvsMC = new TH2F("hnTracksLeadingAODvsMC","",20,0.,20,20,-0.5,19.5);
2042 fhnTracksLeadingAODvsMC->SetXTitle("nTracks MC");
2043 fhnTracksLeadingAODvsMC->SetYTitle("nTracks AOD");
2044 fhnTracksLeadingAODvsMC->Sumw2();
2045 fListOfHistos->Add(fhnTracksLeadingAODvsMC);
2046
2047 fhLeadingdRAODMC = new TH1F("hLeadingdRAODMC","",40,0.,4);
2048 fhLeadingdRAODMC->SetXTitle("#Delta R");
2049 fhLeadingdRAODMC->Sumw2();
2050 fListOfHistos->Add(fhLeadingdRAODMC);
2051
2052
2053 fhLeadingPtAODvsMCdRcut = new TH2F("hLeadingPtAODvsMCdRcut","",40,0,20,40,0.,20.);
2054 fhLeadingPtAODvsMCdRcut->SetXTitle("p_{T} MC (GeV/c)");
2055 fhLeadingPtAODvsMCdRcut->SetYTitle("p_{T} AOD (GeV/c)");
2056 fhLeadingPtAODvsMCdRcut->Sumw2();
2057 fListOfHistos->Add(fhLeadingPtAODvsMCdRcut);
2058
2059
2060 fhdnTracksVsdPtLeadingAODMC = new TH2F("hdnTracksVsdPtLeadingAODMC","",40,-10.,10,20,-10.5,9.5);
2061 fhdnTracksVsdPtLeadingAODMC->SetXTitle("#Delta Pt AOD-MC (GeV/c)");
2062 fhdnTracksVsdPtLeadingAODMC->SetYTitle("#Delta N AOD-MC");
2063 fhdnTracksVsdPtLeadingAODMC->Sumw2();
2064 fListOfHistos->Add(fhdnTracksVsdPtLeadingAODMC);
2065
2066 fhnTracksJetVsPtAOD = new TH2F("hnTracksJetVsPtAOD","",50,0,50,20,-0.5,19.5);
2067 fhnTracksJetVsPtAOD->SetXTitle("p_{T} (GeV/c)");
2068 fhnTracksJetVsPtAOD->SetYTitle("nTracks");
2069 fhnTracksJetVsPtAOD->Sumw2();
2070 fListOfHistos->Add(fhnTracksJetVsPtAOD);
2071
2072 fhnTracksJetVsPtAODquarkEv = new TH2F("hnTracksJetVsPtAODquarkEv","",50,0,50,20,-0.5,19.5);
2073 fhnTracksJetVsPtAODquarkEv->SetXTitle("p_{T} (GeV/c)");
2074 fhnTracksJetVsPtAODquarkEv->SetYTitle("nTracks");
2075 fhnTracksJetVsPtAODquarkEv->Sumw2();
2076 fListOfHistos->Add(fhnTracksJetVsPtAODquarkEv);
2077
2078 fhRadiusJetVsPtAOD = new TH2F("hRadiusJetVsPtAOD","",50,0,50,10,0.,1);
2079 fhRadiusJetVsPtAOD->SetXTitle("p_{T} (GeV/c)");
2080 fhRadiusJetVsPtAOD->SetYTitle("radius");
2081 fhRadiusJetVsPtAOD->Sumw2();
2082 fListOfHistos->Add(fhRadiusJetVsPtAOD);
2083
2084 fhnTracksJetVsPtMC = new TH2F("hnTracksJetVsPtMC","",50,0,50,20,-0.5,19.5);
2085 fhnTracksJetVsPtMC->SetXTitle("p_{T} (GeV/c)");
2086 fhnTracksJetVsPtMC->SetYTitle("nTracks");
2087 fhnTracksJetVsPtMC->Sumw2();
2088 fListOfHistos->Add(fhnTracksJetVsPtMC);
2089
2090 fhnTracksJetVsPtMCquarkEv = new TH2F("hnTracksJetVsPtMCquarkEv","",50,0,50,20,-0.5,19.5);
2091 fhnTracksJetVsPtMCquarkEv->SetXTitle("p_{T} (GeV/c)");
2092 fhnTracksJetVsPtMCquarkEv->SetYTitle("nTracks");
2093 fhnTracksJetVsPtMCquarkEv->Sumw2();
2094 fListOfHistos->Add(fhnTracksJetVsPtMCquarkEv);
2095
2096 fhRadiusJetVsPtMC = new TH2F("hRadiusJetVsPtMC","",50,0,50,10,0.,1);
2097 fhRadiusJetVsPtMC->SetXTitle("p_{T} (GeV/c)");
2098 fhRadiusJetVsPtMC->SetYTitle("radius");
2099 fhRadiusJetVsPtMC->Sumw2();
2100 fListOfHistos->Add(fhRadiusJetVsPtMC);
2101
2102 fhnTracksJetVsPtMCK0 = new TH2F("hnTracksJetVsPtMCK0","",50,0,50,20,-0.5,19.5);
2103 fhnTracksJetVsPtMCK0->SetXTitle("p_{T} (GeV/c)");
2104 fhnTracksJetVsPtMCK0->SetYTitle("nTracks");
2105 fhnTracksJetVsPtMCK0->Sumw2();
2106 fListOfHistos->Add(fhnTracksJetVsPtMCK0);
2107
2108
2109 fhnTracksJetVsPtMCK0quarkEv = new TH2F("hnTracksJetVsPtMCK0quarkEv","",50,0,50,20,-0.5,19.5);
2110 fhnTracksJetVsPtMCK0quarkEv->SetXTitle("p_{T} (GeV/c)");
2111 fhnTracksJetVsPtMCK0quarkEv->SetYTitle("nTracks");
2112 fhnTracksJetVsPtMCK0quarkEv->Sumw2();
2113 fListOfHistos->Add(fhnTracksJetVsPtMCK0quarkEv);
2114
2115 fhRadiusJetVsPtMCK0 = new TH2F("hRadiusJetVsPtMCK0","",50,0,50,10,0.,1);
2116 fhRadiusJetVsPtMCK0->SetXTitle("p_{T} (GeV/c)");
2117 fhRadiusJetVsPtMCK0->SetYTitle("radius");
2118 fhRadiusJetVsPtMCK0->Sumw2();
2119 fListOfHistos->Add(fhRadiusJetVsPtMCK0);
2120
2121 fhnTracksJetVsPtAODK0 = new TH2F("hnTracksJetVsPtAODK0","",50,0,50,20,-0.5,19.5);
2122 fhnTracksJetVsPtAODK0->SetXTitle("p_{T} (GeV/c)");
2123 fhnTracksJetVsPtAODK0->SetYTitle("nTracks AODK0");
2124 fhnTracksJetVsPtAODK0->Sumw2();
2125 fListOfHistos->Add(fhnTracksJetVsPtAODK0);
2126
2127 fhnTracksJetVsPtAODK0quarkEv = new TH2F("hnTracksJetVsPtAODK0quarkEv","",50,0,50,20,-0.5,19.5);
2128 fhnTracksJetVsPtAODK0quarkEv->SetXTitle("p_{T} (GeV/c)");
2129 fhnTracksJetVsPtAODK0quarkEv->SetYTitle("nTracks AODK0quarkEv");
2130 fhnTracksJetVsPtAODK0quarkEv->Sumw2();
2131 fListOfHistos->Add(fhnTracksJetVsPtAODK0quarkEv);
2132
2133 fhRadiusJetVsPtAODK0 = new TH2F("hRadiusJetVsPtAODK0","",50,0,50,10,0.,1);
2134 fhRadiusJetVsPtAODK0->SetXTitle("p_{T} (GeV/c)");
2135 fhRadiusJetVsPtAODK0->SetYTitle("radius");
2136 fhRadiusJetVsPtAODK0->Sumw2();
2137 fListOfHistos->Add(fhRadiusJetVsPtAODK0);
2138
2139 fhnTracksJetVsPtAODpKch = new TH2F("hnTracksJetVsPtAODpKch","",50,0,50,20,-0.5,19.5);
2140 fhnTracksJetVsPtAODpKch->SetXTitle("p_{T} (GeV/c)");
2141 fhnTracksJetVsPtAODpKch->SetYTitle("nTracks AODpKch");
2142 fhnTracksJetVsPtAODpKch->Sumw2();
2143 fListOfHistos->Add(fhnTracksJetVsPtAODpKch);
2144
2145 fhRadiusJetVsPtAODpKch = new TH2F("hRadiusJetVsPtAODpKch","",50,0,50,20,-0.5,19.5);
2146 fhRadiusJetVsPtAODpKch->SetXTitle("p_{T} (GeV/c)");
2147 fhRadiusJetVsPtAODpKch->SetYTitle("Radius AODpKch");
2148 fhRadiusJetVsPtAODpKch->Sumw2();
2149 fListOfHistos->Add(fhRadiusJetVsPtAODpKch);
2150
2151 fhPythiaProcess = CreatePythiaIDhisto("hPythiaProcess");
2152 fListOfHistos->Add(fhPythiaProcess);
2153
2154 fhPythiaProcessK0 = CreatePythiaIDhisto("hPythiaProcessK0");
2155 fListOfHistos->Add(fhPythiaProcessK0);
2156
2157 fhPythiaProcessKch = CreatePythiaIDhisto("hPythiaProcessKch");
2158 fListOfHistos->Add(fhPythiaProcessKch);
2159
2160 fhPythiaProcessp = CreatePythiaIDhisto("hPythiaProcessp");
2161 fListOfHistos->Add(fhPythiaProcessp);
2162
2163 fhPythiaProcesspbar = CreatePythiaIDhisto("hPythiaProcesspbar");
2164 fListOfHistos->Add(fhPythiaProcesspbar);
2165
2166 fh1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
2167 fh1Xsec->SetXTitle("<#sigma>");
2168 fh1Xsec->Sumw2();
2169 fListOfHistos->Add( fh1Xsec );
2170
2171 fh1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
2172 fh1Trials->SetXTitle("#sum{ntrials}");
2173 fh1Trials->Sumw2();
2174 fListOfHistos->Add( fh1Trials );
2175
2176// fSettingsTree = new TTree("JetChemAnalysisSettings","Analysis Settings");
2177// fSettingsTree->Branch("fUseLOConeJets",&fUseLOConeJets,"UseLOConeJets/O");
2178// fSettingsTree->Branch("fUseLOConeMCJets",&fUseLOConeMCJets,"UseLOConeMCJets/O");
2179// fSettingsTree->Branch("fUsePythiaJets",&fUsePythiaJets,"UsePythiaJets/O");
2180// fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
2181// fSettingsTree->Branch("fTrackPtCutJF", &fTrackPtCutJF,"TrackPt/D");
2182// fSettingsTree->Branch("fFilterBitJF",&fFilterBitJF,"FilterBitJF/i");
2183// fSettingsTree->Branch("fRequireITSRefitJF",&fRequireITSRefitJF,"RequireITSRefitJF/O");
2184// fSettingsTree->Branch("fRejectK0TracksJF",&fRejectK0TracksJF,"RejectK0TracksJF/O");
2185// fSettingsTree->Branch("fJetPtCut",&fJetPtCut,"JetPtCut/D");
2186// fSettingsTree->Branch("fJetEtaCut",&fJetEtaCut,"JetEtaCut/D");
2187// fSettingsTree->Branch("fFilterBit",&fFilterBit,"FilterBit/i");
2188// fSettingsTree->Branch("fTrackPtCut",&fTrackPtCut,"TrackPtCut/D");
2189// fSettingsTree->Branch("fTrackEtaCut",&fTrackEtaCut,"TrackEtaCut/D");
2190// fSettingsTree->Branch("fUseOnFlyV0s",&fUseOnFlyV0s,"UseOnFlyV0s/O");
2191// fSettingsTree->Branch("fCutnSigdEdx",&fCutnSigdEdx,"CutnSigdEdx/D");
2192// fSettingsTree->Branch("fUseAODMCTracksForUE",&fUseAODMCTracksForUE,"UseAODMCTracksForUE/O");
2193// fSettingsTree->Branch("fAreaReg",&fAreaReg,"AreaReg/D");
2194// fSettingsTree->Branch("fAvgTrials",&fAvgTrials,"AvgTrials/D");
2195 // fListOfHistos->Add(fSettingsTree);
2196
2197}
2198
2199//____________________________________________________________________
2200void AliAnalysisTaskJetChem::Terminate(Option_t */*option*/){
2201
2202 // Terminate analysis
2203 //
2204
2205 WriteSettings();
2206
2207 // Update pointers reading them from the output slot
2208 fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
2209 if( !fListOfHistos ) {
2210 AliError("Histogram List is not available");
2211 return;
2212 }
2213
2214 fhLeadingPt = (TH1F*)fListOfHistos->FindObject("hLeadingPt");
2215 fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->FindObject("hRegionSumPtMaxVsEt");
2216 fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->FindObject("hRegionSumPtMinVsEt");
2217 fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->FindObject("hRegionMultMaxVsEt");
2218 fhRegionMultMinVsEt = (TH1F*)fListOfHistos->FindObject("hRegionMultMinVsEt");
2219
2220
2221 fhRegionSumPtMaxVsEt->Divide(fhLeadingPt);
2222 fhRegionSumPtMinVsEt->Divide(fhLeadingPt);
2223 fhRegionMultMaxVsEt->Divide(fhLeadingPt);
2224 fhRegionMultMinVsEt->Divide(fhLeadingPt);
2225
2226
2227 //Get Normalization
2228 fh1Xsec = (TProfile*) fListOfHistos->FindObject("fh1Xsec");
2229 fh1Trials = (TH1F*)fListOfHistos->FindObject("fh1Trials");
2230
2231 //std::cout<<" fh1Xsec "<<fh1Xsec<<" fh1Trials "<<fh1Trials<<std::endl;
2232
2233 if(fh1Xsec && fh1Trials){
2234
2235 Double_t xsec = fh1Xsec->GetBinContent(1);
2236 Double_t ntrials = fh1Trials->GetBinContent(1);
2237 Double_t normFactor = xsec/ntrials;
2238 Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
2239 fhLeadingPt->Scale(normFactor);
2240 }
2241
2242
2243 if(fDebug > 1) AliInfo("End analysis");
2244}
2245
2246// ----------------------------------------------------------------------------
2247
2248void AliAnalysisTaskJetChem::WriteSettings(){
2249
2250 // fill settings tree
2251 // not on GRID !
2252
2253 //fSettingsTree->Fill();
2254
2255}
2256
2257// ---------------------------------------------------------------------------
2258
2259Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const {
2260
2261 // K0 mass ? Use STAR params for bin counting & mass fit
2262
2263 const Double_t massK0 = 0.497; // from fits
2264 const Double_t sigmaK0 = 0.0046;
2265 const Double_t nSigmaSignal = 3.5; // STAR parameters for bin counting
2266
2267 if((massK0-nSigmaSignal*sigmaK0)<=mass &&
2268 mass<=(massK0 + nSigmaSignal*sigmaK0)) return kTRUE;
2269
2270 return kFALSE;
2271}
2272
2273// ---------------------------------------------------------------------------
2274
2275Bool_t AliAnalysisTaskJetChem::IsLambdaInvMass(const Double_t mass) const{
2276
2277 // Lambda mass ?
2278
2279
2280 if(1.1<mass && mass<1.13) return kTRUE; // FIXME - adjust range from fit
2281 return kFALSE;
2282}
2283
2284// ---------------------------------------------------------------------------
2285
2286Bool_t AliAnalysisTaskJetChem::IsAcceptedDCAK0(const Double_t dca) const{
2287
2288 // DCA cut
2289
2290 return kTRUE; // FIXME - adjust cut
2291}
2292
2293
2294// ---------------------------------------------------------------------------
2295
2296Bool_t AliAnalysisTaskJetChem::IsAcceptedDCALambda(const Double_t dca) const {
2297
2298 // DCA cut
2299
2300 return kTRUE; // FIXME - adjust cut
2301}
2302
2303// ---------------------------------------------------------------------------
2304
2305Bool_t AliAnalysisTaskJetChem::IsAccepteddEdx(const Double_t mom,const Double_t signal, AliPID::EParticleType n, const Double_t cutnSig) const{
2306
2307 // apply TPC dE/dx cut similar as in AliTPCpidESD
2308 // note: AliTPCidESD uses inner track param for momentum - not avaiable on AOD,
2309 // so we use global track momentum
2310 // should use separate parametrisation for MC and data, but probably ALEPH param & 7% resolution used here anyway not the last word
2311
2312
2313 const Double_t kBBMIP(50.);
2314 const Double_t kBBRes(0.07);
2315 //const Double_t kBBRange(5.);
2316 const Double_t kBBp1(0.76176e-1);
2317 const Double_t kBBp2(10.632);
2318 const Double_t kBBp3(0.13279e-4);
2319 const Double_t kBBp4(1.8631);
2320 const Double_t kBBp5(1.9479);
2321
2322 Double_t mass=AliPID::ParticleMass(n);
2323 Double_t betaGamma = mom/mass;
2324
2325 const Float_t kmeanCorrection =0.1;
2326 Double_t bb = AliExternalTrackParam::BetheBlochAleph(betaGamma,kBBp1,kBBp2,kBBp3,kBBp4,kBBp5);
2327 Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
2328 Double_t bethe = bb * meanCorrection; // expected
2329 Double_t sigma = bethe * kBBRes;
2330
2331
2332 Double_t dedx = signal/kBBMIP; // measured
2333
2334 Double_t nSig = (TMath::Abs(dedx - bethe))/sigma;
2335
2336 if(nSig > cutnSig) return kFALSE;
2337
2338 return kTRUE;
2339}
2340
2341// ----------------------------------------------------------------------------
2342
2343void AliAnalysisTaskJetChem::CheckV0s(AliAODJet* jetVect, Int_t maxPtRegionIndex,Bool_t& foundK0){
2344
2345 // loop over AOD V0s, fill masses etc
2346
2347 Int_t nV0 = fAOD->GetNumberOfV0s();
2348 fhNV0s->Fill(nV0);
2349
2350 for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
2351
2352 AliAODv0* v0 = fAOD->GetV0(i);
2353
2354 Bool_t isOnFly = v0->GetOnFlyStatus();
2355 isOnFly ? fhV0onFly->Fill(1) : fhV0onFly->Fill(0);
2356
2357 if((fUseOnFlyV0s && isOnFly) || (!fUseOnFlyV0s && !isOnFly) ){ // 'offline' V0s : using vertex tracks, on-fly: take decay vertex into account during mom fit
2358
2359 Double_t massK0 = v0->MassK0Short();
2360 Double_t massLambda = v0->MassLambda();
2361 Double_t massAntiLambda = v0->MassAntiLambda();
2362 Double_t radiusV0 = v0->RadiusV0();
2363 Double_t etaV0 = v0->Eta();
2364 Double_t fDCAV0ToVertex = v0->DcaV0ToPrimVertex();
2365 Double_t fDCADaughters = v0->DcaV0Daughters();
2366
2367
2368 Double_t v0Mom[3];
2369 v0->PxPyPz(v0Mom);
2370 TVector3 v0MomVect(v0Mom);
2371
2372 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2373 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2374
2375 Double_t posMom[3], negMom[3];
2376 trackPos->PxPyPz(posMom);
2377 trackNeg->PxPyPz(negMom);
2378
2379 TVector3 posMomVec(posMom);
2380 TVector3 negMomVec(negMom);
2381
2382 Double_t dROpanV0 = posMomVec.DeltaR(negMomVec);
2383
2384 AliAODPid* aodPidPos = trackPos->GetDetPid();
2385 AliAODPid* aodPidNeg = trackNeg->GetDetPid();
2386
2387 Double_t dEdxPos = aodPidPos->GetTPCsignal();
2388 Double_t dEdxNeg = aodPidNeg->GetTPCsignal();
2389
2390 Double_t momPos = trackPos->P();
2391 Double_t momNeg = trackNeg->P();
2392
2393 fhdEdxVsMomV0->Fill(momPos,dEdxPos);
2394 fhdEdxVsMomV0->Fill(momNeg,dEdxNeg);
2395
2396 Double_t pV0 = TMath::Sqrt(v0->Ptot2V0());
2397 Double_t ptV0 = TMath::Sqrt(v0->Pt2V0());
2398
2399 fhV0InvMassK0->Fill(massK0);
2400 fhV0InvMassLambda->Fill(massLambda);
2401 fhV0InvMassAntiLambda->Fill(massAntiLambda);
2402
2403 fhV0DCADaughters->Fill(fDCADaughters);
2404 fhV0Radius->Fill(radiusV0);
2405 fhV0DCAToVertex->Fill(fDCAV0ToVertex);
2406 fhdNdptV0->Fill(ptV0);
2407
2408
2409 // K0 signal before cuts
2410
2411 if(IsK0InvMass(massK0)){
2412
2413 fhdNdptK0->Fill(ptV0);
2414
2415 fhPtVsEtaK0->Fill(etaV0,ptV0);
2416
2417 Double_t dRV0MC = AssociateV0MC(&v0MomVect,310); // K0
2418 if(dRV0MC < 0) dRV0MC = 0.99;
2419 fhdRV0MC->Fill(dRV0MC);
2420 }
2421
2422 if(IsAcceptedDCAK0(fDCAV0ToVertex)){
2423
2424 fhV0InvMassK0DCA->Fill(massK0);
2425 if(IsK0InvMass(massK0)) fhdNdptK0DCA->Fill(ptV0);
2426
2427
2428 if(IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)) fhdEdxVsMomV0pidEdx->Fill(momPos,dEdxPos);
2429 if(IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) fhdEdxVsMomV0pidEdx->Fill(-1*momNeg,dEdxNeg);
2430
2431 if(IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) &&
2432 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2433
2434 fhV0InvMassK0DCAdEdx->Fill(massK0);
2435
2436 if(IsK0InvMass(massK0)){
2437 fhdNdptK0DCAdEdx->Fill(ptV0);
2438 fhV0DCAToVertexK0->Fill(fDCAV0ToVertex);
2439 fhdROpanK0VsPt->Fill(ptV0,dROpanV0);
2440 foundK0 = kTRUE;
2441 }
2442 }
2443
2444
2445 // AOD pid - cuts strongly into signal
2446
2447 AliAODTrack::AODTrkPID_t mpPIDNeg = trackNeg->GetMostProbablePID();
2448 AliAODTrack::AODTrkPID_t mpPIDPos = trackPos->GetMostProbablePID();
2449
2450 if(mpPIDNeg == AliAODTrack::kPion) fhdEdxVsMomV0piPID->Fill(momPos,dEdxPos);
2451 if(mpPIDPos == AliAODTrack::kPion) fhdEdxVsMomV0piPID->Fill(-1*momNeg,dEdxNeg);
2452
2453
2454 if( (mpPIDNeg == AliAODTrack::kPion) && (mpPIDPos == AliAODTrack::kPion) ){
2455
2456 fhV0InvMassK0DCAPID->Fill(massK0);
2457 }
2458 }
2459
2460
2461 if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
2462
2463 fhV0InvMassK0Lambda->Fill(massK0);
2464 }
2465
2466
2467 // Lambda
2468
2469 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2470 IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2471 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2472
2473 fhV0InvMassLambdaDCAdEdx->Fill(massLambda);
2474 }
2475
2476 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2477 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2478 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
2479
2480 fhV0InvMassAntiLambdaDCAdEdx->Fill(massAntiLambda);
2481 }
2482
2483
2484 // jet events
2485
2486 if(jetVect){
2487
2488 Int_t regionV0vect = IsTrackInsideRegion(jetVect,&v0MomVect);
2489
2490 // calc xi
2491 Double_t jetpt = jetVect->Pt();
2492 Double_t dPhiJetV0 = (jetVect->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
2493 Double_t dRJetV0 = (jetVect->MomentumVector()->Vect()).DeltaR(v0MomVect);
2494
2495 Double_t z = pV0/jetpt;
2496 Double_t xi = TMath::Log(1/z);
2497
2498 fhdPhiJetV0->Fill(dPhiJetV0);
2499
2500 // K0
2501
2502 if(IsAcceptedDCAK0(fDCAV0ToVertex) &&
2503 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) &&
2504 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2505
2506 fhV0InvMassK0JetEvt->Fill(massK0);
2507
2508 if(IsK0InvMass(massK0)){
2509 fhdNdptK0JetEvt->Fill(ptV0);
2510 fhdNdzK0->Fill(z);
2511 fhdNdxiK0->Fill(xi);
2512
2513 fhdPhiJetK0->Fill(dPhiJetV0);
2514 fhdRJetK0->Fill(dRJetV0);
2515 }
2516 }
2517
2518
2519 // Lambda
2520
2521 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2522 IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2523 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2524
2525 fhV0InvMassLambdaJetEvt->Fill(massLambda);
2526 }
2527
2528 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2529 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2530 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
2531
2532 fhV0InvMassAntiLambdaJetEvt->Fill(massAntiLambda);
2533 }
2534
2535 // fill histos max region
2536
2537 if(regionV0vect != 0 && regionV0vect == maxPtRegionIndex){ // max region
2538
2539 // K0
2540
2541 if(IsAcceptedDCAK0(fDCAV0ToVertex) &&
2542 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) &&
2543 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){ // K0 cuts
2544
2545 fhV0InvMassK0Max->Fill(massK0);
2546
2547 if(IsK0InvMass(massK0)){
2548 fhdNdzK0Max->Fill(z);
2549 fhdNdxiK0Max->Fill(xi);
2550 fhdNdptK0Max->Fill(ptV0);
2551 }
2552 }
2553
2554 // Lambda
2555
2556 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2557 IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2558 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2559
2560 fhV0InvMassLambdaMax->Fill(massLambda);
2561 }
2562
2563 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2564 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2565 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
2566
2567 fhV0InvMassAntiLambdaMax->Fill(massAntiLambda);
2568 }
2569
2570 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2571 (IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2572 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) ||
2573 (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2574 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx))){
2575
2576 if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
2577 fhdNdzLambdaMax->Fill(z);
2578 fhdNdxiLambdaMax->Fill(xi);
2579 fhdNdptLambdaMax->Fill(ptV0);
2580 }
2581 }
2582 }
2583
2584 // fill histos min region
2585
2586 if(regionV0vect != 0 && regionV0vect != maxPtRegionIndex){ // min region
2587
2588 // K0
2589
2590 if(IsAcceptedDCAK0(fDCAV0ToVertex) &&
2591 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) &&
2592 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2593
2594 fhV0InvMassK0Min->Fill(massK0);
2595
2596 if(IsK0InvMass(massK0)){
2597 fhdNdzK0Min->Fill(z);
2598 fhdNdxiK0Min->Fill(xi);
2599 fhdNdptK0Min->Fill(ptV0);
2600 }
2601 }
2602
2603
2604 // Lambda
2605
2606
2607 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2608 IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2609 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2610
2611 fhV0InvMassLambdaMin->Fill(massLambda);
2612 }
2613
2614 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2615 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2616 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
2617
2618 fhV0InvMassAntiLambdaMin->Fill(massAntiLambda);
2619 }
2620
2621 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2622 (IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2623 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) ||
2624 (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2625 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx))){
2626
2627 if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
2628 fhdNdzLambdaMin->Fill(z);
2629 fhdNdxiLambdaMin->Fill(xi);
2630 fhdNdptLambdaMin->Fill(ptV0);
2631 }
2632 }
2633 }
2634
2635
2636 // jet region
2637
2638 if(regionV0vect == 0){ // jet region
2639
2640 //Double_t dRJetV0 = jetVect->DeltaR(v0MomVect);
2641
2642 if(dRJetV0 <= fConeRadius){
2643
2644 // fill histos jet region
2645
2646 // K0
2647
2648 if(IsAcceptedDCAK0(fDCAV0ToVertex) &&
2649 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) &&
2650 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){ // K0 cuts
2651
2652 fhV0InvMassK0Jet->Fill(massK0);
2653
2654 if(IsK0InvMass(massK0)){
2655 fhdNdzK0Jet->Fill(z);
2656 fhdNdxiK0Jet->Fill(xi);
2657 fhdNdptK0Jet->Fill(ptV0);
2658 }
2659 }
2660
2661
2662 // Lambda
2663
2664 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2665 IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2666 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
2667
2668 fhV0InvMassLambdaJet->Fill(massLambda);
2669 }
2670 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2671 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2672 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
2673
2674 fhV0InvMassAntiLambdaJet->Fill(massAntiLambda);
2675 }
2676
2677 if(IsAcceptedDCALambda(fDCAV0ToVertex) &&
2678 (IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) &&
2679 IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) ||
2680 (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) &&
2681 IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx))){
2682
2683 if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
2684 fhdNdzLambdaJet->Fill(z);
2685 fhdNdxiLambdaJet->Fill(xi);
2686 fhdNdptLambdaJet->Fill(ptV0);
2687 }
2688 }
2689 }
2690 }
2691 }
2692 }
2693 }
2694}
2695
2696// ----------------------------------------------------------------------------
2697
2698void AliAnalysisTaskJetChem::CheckMCParticles(AliAODJet* jetVect,Int_t maxPtRegionIndex,Bool_t& isK0Event){
2699
2700 // histos for generated MC
2701
2702 TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
2703
2704 if(!farray){
2705 AliInfo("no mcparticles branch");
2706 return;
2707 }
2708
2709 Int_t ntrks = farray->GetEntries();
2710 if (fDebug>1) AliInfo(Form("check MC particles, tracks %d \n",ntrks));
2711
2712 Int_t pythiaPID = GetPythiaProcessID();
2713
2714 isK0Event = kFALSE;
2715 Bool_t ispEvent = kFALSE;
2716 Bool_t ispBarEvent = kFALSE;
2717 Bool_t isKchEvent = kFALSE;
2718
2719 for(Int_t i =0 ; i<ntrks; i++){ // mc tracks loop
2720
2721 AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
2722
2723 //Cuts
2724 if (!(mctrk->IsPhysicalPrimary())) continue;
2725
2726 if ((mctrk->Pt() < fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut )) continue;
2727
2728 // OB PID histo
2729 Int_t pdg = mctrk->GetPdgCode();
2730
2731 FillPIDhisto(fhPIDMC,pdg);
2732
2733 if(pdg == 111) fhdNdptpi0MC->Fill(mctrk->Pt());
2734 if(pdg == 22) fhdNdptgammaMC->Fill(mctrk->Pt());
2735 if(TMath::Abs(pdg) == 211) fhdNdptchPiMC->Fill(mctrk->Pt());
2736 if(pdg == 310) fhdNdptK0MC->Fill(mctrk->Pt());
2737 if(TMath::Abs(pdg) == 321) fhdNdptchKMC->Fill(mctrk->Pt());
2738 if(pdg == 2212) fhdNdptpMC->Fill(mctrk->Pt());
2739 if(pdg == -2212) fhdNdptpBarMC->Fill(mctrk->Pt());
2740 if(TMath::Abs(pdg) == 3122) fhdNdptLambdaLambdaBarMC->Fill(mctrk->Pt());
2741 if(pdg == 3332) fhdNdptOmegaMC->Fill(mctrk->Pt());
2742 if(pdg == -3332) fhdNdptOmegaBarMC->Fill(mctrk->Pt());
2743
2744 if(pdg == 310) isK0Event = kTRUE;
2745 if(TMath::Abs(pdg) == 321) isKchEvent = kTRUE;
2746 if(pdg == 2212) ispEvent = kTRUE;
2747 if(pdg == -2212) ispBarEvent = kTRUE;
2748
2749
2750 if(jetVect){ // jet event
2751
2752 FillPIDhisto(fhPIDMCAll,pdg);
2753
2754 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
2755
2756 Int_t region = IsTrackInsideRegion(jetVect, &partVect );
2757
2758 //Double_t deltaPhi = jetVect->DeltaPhi(partVect); //+k270rad;
2759 Double_t deltaPhi = (jetVect->MomentumVector()->Vect()).DeltaPhi(partVect);
2760 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
2761 //Double_t deltaR = jetVect->DeltaR(partVect); //+k270rad;
2762 Double_t deltaR = (jetVect->MomentumVector()->Vect()).DeltaR(partVect);
2763
2764
2765 // calc xi
2766 Double_t jetpt = jetVect->Pt();
2767 Double_t pV0 = partVect.Mag();
2768 Double_t ptV0 = partVect.Pt();
2769
2770 Double_t z = pV0 / jetpt;
2771 Double_t xi = TMath::Log(1/z);
2772
2773 if(!(mctrk->Charge() == 0 || mctrk->Charge()==-99)) fhdNdxiMC->Fill(xi);
2774
2775 if(TMath::Abs(pdg) == 310){ // K0
2776 fhdPhiJetK0MC->Fill(deltaPhi);
2777 fhdRJetK0MC->Fill(deltaR);
2778
2779 fhdNdxiK0MC->Fill(xi);
2780 fhdNdzK0MC->Fill(z);
2781 fhdNdptK0MCJetEvt->Fill(ptV0);
2782 }
2783
2784
2785 if(region != 0 && region == maxPtRegionIndex){ // max region
2786
2787 if(TMath::Abs(pdg) == 310){ // K0
2788
2789 fhdNdptK0MCMax->Fill(ptV0);
2790 }
2791 if(TMath::Abs(pdg) == 3122){ // lambda / anti-lambda
2792
2793 fhdNdptLambdaLambdaBarMCMax->Fill(ptV0);
2794 }
2795
2796 }
2797
2798 if(region != 0 && region != maxPtRegionIndex){ // min region
2799
2800 FillPIDhisto(fhPIDMCMin,pdg);
2801
2802 if(TMath::Abs(pdg) == 310){ // K0
2803
2804 fhdNdptK0MCMin->Fill(ptV0);
2805 }
2806 if(TMath::Abs(pdg) == 3122){ // lambda / anti-lambda
2807
2808 fhdNdptLambdaLambdaBarMCMin->Fill(ptV0);
2809 }
2810 if(TMath::Abs(pdg) == 211){ // ch pi
2811
2812 fhdNdptchPiMCMin->Fill(ptV0);
2813 }
2814 if(TMath::Abs(pdg) == 321){ // ch K
2815
2816 fhdNdptchKMCMin->Fill(ptV0);
2817 }
2818 if(pdg == 2212){ // p
2819
2820 fhdNdptpMCMin->Fill(ptV0);
2821 }
2822 if(pdg == -2212){ // p
2823
2824 fhdNdptpBarMCMin->Fill(ptV0);
2825 }
2826
2827 if(pdg == 3332) fhdNdptOmegaMCMin->Fill(mctrk->Pt());
2828 if(pdg == -3332) fhdNdptOmegaBarMCMin->Fill(mctrk->Pt());
2829 }
2830
2831 if(region == 0){ // jet region ?
2832
2833 FillPIDhisto(fhPIDMCJet,pdg);
2834
2835 //Double_t dRJetV0 = jetVect->DeltaR(partVect);
2836 Double_t dRJetV0 = (jetVect->MomentumVector()->Vect()).DeltaR(partVect);
2837
2838 if(dRJetV0 <= fConeRadius){
2839
2840 if(TMath::Abs(pdg) == 310){ // K0
2841
2842 fhdNdptK0MCJet->Fill(ptV0);
2843 fhdNdxiK0MCJet->Fill(xi);
2844 fhdNdzK0MCJet->Fill(z);
2845 }
2846 if(TMath::Abs(pdg) == 3122){ // lambda / anti-lambda
2847
2848 fhdNdptLambdaLambdaBarMCJet->Fill(ptV0);
2849 }
2850 if(TMath::Abs(pdg) == 211){ // ch pi
2851
2852 fhdNdptchPiMCJet->Fill(ptV0);
2853 }
2854 if(TMath::Abs(pdg) == 321){ // ch K
2855
2856 fhdNdptchKMCJet->Fill(ptV0);
2857 }
2858 if(pdg == 2212){ // p/pbar
2859
2860 fhdNdptpMCJet->Fill(ptV0);
2861 }
2862 if(pdg == -2212){ // p/pbar
2863
2864 fhdNdptpBarMCJet->Fill(ptV0);
2865 }
2866 }
2867 }
2868 }
2869 }
2870
2871
2872 FillPythiaIDhisto(fhPythiaProcess,pythiaPID);
2873 if(isK0Event) FillPythiaIDhisto(fhPythiaProcessK0,pythiaPID);
2874 if(isKchEvent) FillPythiaIDhisto(fhPythiaProcessKch,pythiaPID);
2875 if(ispEvent) FillPythiaIDhisto(fhPythiaProcessp,pythiaPID);
2876 if(ispBarEvent) FillPythiaIDhisto(fhPythiaProcesspbar,pythiaPID);
2877
2878}
2879
2880// ----------------------------------------------------------------------------
2881
2882Double_t AliAnalysisTaskJetChem::AssociateV0MC(const TVector3* V0Mom,const Int_t pdg){
2883
2884 // find closest MC gen. particle for V0 vector
2885
2886 TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
2887
2888 if(!farray){
2889 AliInfo("no mcparticles branch");
2890 return -1;
2891 }
2892
2893 Double_t dRmin = -1;
2894
2895 Int_t ntrks = farray->GetEntries();
2896
2897 for(Int_t i =0 ; i<ntrks; i++){
2898
2899 AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
2900
2901 //Cuts
2902 if (!(mctrk->IsPhysicalPrimary())) continue;
2903
2904 Int_t pdgtrk = mctrk->GetPdgCode();
2905
2906 if(pdgtrk != pdg) continue;
2907
2908 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
2909
2910 Double_t dR = V0Mom->DeltaR(partVect);
2911
2912 if(dRmin<0) dRmin = dR; // initialize
2913
2914 if(dR < dRmin) dRmin = dR;
2915
2916 }
2917
2918 return dRmin;
2919}
2920
2921// ----------------------------------------------------------------------------
2922
2923void AliAnalysisTaskJetChem::CompLeadingJets(AliAODJet* jetLeadingAOD,AliAODJet* jetLeadingMC,const Int_t pythiaPID,
2924 const Bool_t foundK0AOD,const Bool_t foundK0MC){
2925
2926 // leading jet properties
2927
2928 Double_t ptLeadingAOD = -1;
2929 Double_t etaLeadingAOD = -1;
2930 Double_t phiLeadingAOD = -1;
2931 Int_t nTracksLeadingAOD = -1;
2932
2933 Double_t ptLeadingMC = -1;
2934 Double_t etaLeadingMC = -1;
2935 Double_t phiLeadingMC = -1;
2936 Int_t nTracksLeadingMC = -1;
2937
2938 if(jetLeadingAOD){
2939 if(jetLeadingAOD->Pt()>fJetPtCut && TMath::Abs(jetLeadingAOD->Eta())<fJetEtaCut){
2940
2941 ptLeadingAOD = jetLeadingAOD->Pt();
2942 etaLeadingAOD = jetLeadingAOD->Eta();
2943 phiLeadingAOD = jetLeadingAOD->Phi();
2944 nTracksLeadingAOD = jetLeadingAOD->GetRefTracks()->GetEntriesFast();
2945
2946 Double_t radiusAOD = GetJetRadius(jetLeadingAOD,0.8);
2947 fhnTracksJetVsPtAOD->Fill(ptLeadingAOD,nTracksLeadingAOD);
2948 fhRadiusJetVsPtAOD->Fill(ptLeadingAOD,radiusAOD);
2949 if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtAODquarkEv->Fill(ptLeadingAOD,nTracksLeadingAOD);
2950
2951 if(foundK0AOD){
2952
2953 fhnTracksJetVsPtAODK0->Fill(ptLeadingAOD,nTracksLeadingAOD);
2954 if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtAODK0quarkEv->Fill(ptLeadingAOD,nTracksLeadingAOD);
2955 fhRadiusJetVsPtAODK0->Fill(ptLeadingAOD,radiusAOD);
2956 }
2957
2958
2959 // check if p/Kch in jet
2960
2961 Bool_t foundpKch = kFALSE;
2962 Int_t nTracksJet = jetLeadingAOD->GetRefTracks()->GetEntriesFast();
2963
2964 for(int i=0; i<nTracksJet; i++){
2965
2966 AliAODTrack* track = (AliAODTrack*) jetLeadingAOD->GetRefTracks()->At(i);
2967
2968 Double_t mom = track->P();
2969
2970 AliAODPid* aodPid = track->GetDetPid();
2971 Double_t dEdx = aodPid->GetTPCsignal();
2972
2973 if(IsAccepteddEdx(mom,dEdx,AliPID::kKaon,fCutnSigdEdx) ||
2974 IsAccepteddEdx(mom,dEdx,AliPID::kProton,fCutnSigdEdx)){
2975
2976 foundpKch = kTRUE;
2977 }
2978 } // track loop
2979
2980 if(foundpKch){
2981 fhnTracksJetVsPtAODpKch->Fill(ptLeadingAOD,nTracksLeadingAOD);
2982 fhRadiusJetVsPtAODpKch->Fill(ptLeadingAOD,radiusAOD);
2983 }
2984 }
2985 }
2986
2987
2988 if(jetLeadingMC){
2989 if(jetLeadingMC->Pt()>fJetPtCut && TMath::Abs(jetLeadingMC->Eta())<fJetEtaCut){
2990
2991 ptLeadingMC = jetLeadingMC->Pt();
2992 etaLeadingMC = jetLeadingMC->Eta();
2993 phiLeadingMC = jetLeadingMC->Phi();
2994 nTracksLeadingMC = jetLeadingMC->GetRefTracks()->GetEntriesFast();
2995
2996 Double_t radiusMC = GetJetRadius(jetLeadingMC,0.8);
2997 fhnTracksJetVsPtMC->Fill(ptLeadingMC,nTracksLeadingMC);
2998 fhRadiusJetVsPtMC->Fill(ptLeadingMC,radiusMC);
2999 if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtMCquarkEv->Fill(ptLeadingMC,nTracksLeadingMC);
3000
3001 if(foundK0MC){
3002
3003 fhnTracksJetVsPtMCK0->Fill(ptLeadingMC,nTracksLeadingMC);
3004 if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtMCK0quarkEv->Fill(ptLeadingMC,nTracksLeadingMC);
3005 fhRadiusJetVsPtMCK0->Fill(ptLeadingMC,radiusMC);
3006 }
3007 }
3008 }
3009
3010 if(jetLeadingAOD && jetLeadingMC){
3011
3012 //std::cout<<" comp: leading jetPt AOD "<<ptLeadingAOD<<" MC "<<ptLeadingMC<<std::endl;
3013 //if(jetLeadingAOD && jetLeadingMC)
3014 //std::cout<<" leading jet eta AOD "<<jetLeadingAOD->Eta()<<" MC "<<jetLeadingMC->Eta()<<std::endl;
3015
3016
3017 if(jetLeadingMC->Pt()>fJetPtCut && jetLeadingAOD->Pt()>fJetPtCut &&
3018 TMath::Abs(jetLeadingMC->Eta())<fJetEtaCut && TMath::Abs(jetLeadingAOD->Eta())<fJetEtaCut){
3019
3020 fhLeadingPtAODvsMC->Fill(ptLeadingMC,ptLeadingAOD);
3021 fhLeadingEtaAODvsMC->Fill(etaLeadingMC,etaLeadingAOD);
3022 fhLeadingPhiAODvsMC->Fill(phiLeadingMC,phiLeadingAOD);
3023 fhnTracksLeadingAODvsMC->Fill(nTracksLeadingMC,nTracksLeadingAOD);
3024
3025 TLorentzVector *mom4MC = jetLeadingMC->MomentumVector();
3026 TLorentzVector *mom4AOD = jetLeadingAOD->MomentumVector();
3027
3028 Double_t dR = mom4MC->DeltaR(*mom4AOD);
3029 fhLeadingdRAODMC->Fill(dR);
3030
3031 if(dR<0.4) fhLeadingPtAODvsMCdRcut->Fill(ptLeadingMC,ptLeadingAOD);
3032
3033 Double_t dPt = ptLeadingAOD - ptLeadingMC;
3034 Double_t dnTracks = nTracksLeadingAOD - nTracksLeadingMC;
3035
3036 fhdnTracksVsdPtLeadingAODMC->Fill(dPt,dnTracks);
3037 }
3038 }
3039
3040}
3041
3042// ---------------------------------------------------------------------------
3043
3044// void AliAnalysisTaskJetChem::CheckK0(){
3045
3046// TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
3047
3048// if(!farray){
3049// AliInfo("no mcparticles branch");
3050// return;
3051// }
3052
3053
3054// Int_t ntrks = farray->GetEntries();
3055
3056// for(Int_t i=0; i<ntrks; i++){ // trk loop
3057
3058// AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
3059
3060// Int_t pdg = mctrk->GetPdgCode();
3061
3062// if(pdg != 310) continue; // K0
3063
3064// cout<<" MC K0: physPrim "<<mctrk->IsPhysicalPrimary()<<endl;
3065
3066// Int_t labelMother = mctrk->GetMother();
3067
3068// cout<<" found K00, label mother "<<labelMother<<endl;
3069
3070// AliAODMCParticle* mctrkMother = NULL;
3071// Int_t pdgMother = -1;
3072
3073// for(Int_t k=0 ; k<ntrks; k++){
3074
3075// mctrkMother = (AliAODMCParticle*)farray->At(k);
3076// if(mctrkMother->GetLabel() == labelMother) break;
3077// }
3078
3079// pdgMother = mctrkMother->GetPdgCode();
3080// cout<<" K0 mother pdg "<<pdgMother<<" GID "<<fpdgdb->ConvertPdgToGeant3(pdgMother)<<" isPrimary "<<mctrkMother->IsPrimary()<<endl;
3081// //cout<<" mother name "<<mctrkMother->GetName()<<endl;
3082
3083// }
3084// }
3085
3086// ---------------------------------------------------------------------------
3087
3088// void AliAnalysisTaskJetChem::CheckK0Stack(){
3089
3090// AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
3091// if (!mcHandler) {
3092// Printf("ERROR: Could not retrieve MC event handler");
3093// return;
3094// }
3095
3096// AliMCEvent* mcEvent = mcHandler->MCEvent();
3097// if (!mcEvent) {
3098// Printf("ERROR: Could not retrieve MC event");
3099// return;
3100// }
3101
3102// AliStack* mcStack = mcEvent->Stack();//Load Stack
3103// AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
3104
3105// Int_t nTracksMC = mcStack->GetNtrack();
3106
3107// Bool_t foundK0 = kFALSE;
3108
3109// for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
3110// //Cuts
3111// if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
3112
3113// TParticle* mctrk = mcStack->Particle(iTracks);
3114
3115// Int_t pdg = mctrk->GetPdgCode();
3116
3117// if ((mctrk->Pt() < 2*fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut )) continue;
3118
3119// if(pdg == 310) foundK0 = kTRUE;
3120
3121// if(pdg != 310) continue; // K0 short
3122
3123// cout<<" check K0 "<<endl;
3124
3125// Int_t indexMother = -999;
3126// TParticle* mctrkMother = mctrk;
3127// TParticle* mctrkFirstGenMother = NULL; // last mother which is not the primary proton
3128
3129// Int_t nAncestors = 0;
3130
3131// while(indexMother != -1){
3132
3133// indexMother = mctrkMother->GetFirstMother();
3134
3135// if(indexMother != -1){
3136// mctrkFirstGenMother = mctrkMother;
3137// mctrkMother = mcStack->Particle(indexMother);
3138// nAncestors++;
3139// }
3140
3141// cout<<" nAncestors "<<nAncestors<<" pdg mother "<<mctrkMother->GetPdgCode()<<" name "<<mctrkMother->GetName()<<endl;
3142// }
3143// cout<<" pdg firstGenMother "<<mctrkFirstGenMother->GetPdgCode()<<" name "<<mctrkFirstGenMother->GetName()<<endl;
3144
3145
3146
3147// cout<<" pythiaGenHeader "<<pythiaGenHeader<<endl;
3148// cout<<" Pythia Process type "<<pythiaGenHeader->ProcessType()<<" ptHard "<<pythiaGenHeader->GetPtHard()<<endl;
3149
3150
3151// fhPythiaProcess->Fill(pythiaGenHeader->ProcessType());
3152// if(foundK0) fhPythiaProcess_K0->Fill(pythiaGenHeader->ProcessType());
3153
3154// //Int_t indexGrandMother = mctrkMother->GetFirstMother();
3155// //cout<<" indexGrandMother "<<indexGrandMother<<endl;
3156
3157// //if(indexGrandMother>-1){
3158// // TParticle* mctrkGrandMother = mcStack->Particle(indexGrandMother);
3159// // cout<<" pdg grandMother "<<mctrkGrandMother->GetPdgCode()<<" name "<<mctrkGrandMother->GetName()<<endl;
3160// // }
3161
3162// }
3163// }
3164
3165// ---------------------------------------------------------------------------------
3166
3167Bool_t AliAnalysisTaskJetChem::IsQuarkHardScatteringEvent(const Int_t PID){
3168
3169 // Pythia Manual sec. 8.2.1 :
3170 // if Pythia PID = 92,93,94 event is diffractive
3171 // if Pythia PID = 13, 28, 68 hard scattering products are gg, qg
3172 // if Pythia PID = 11, 12, 53 hard scattering products are qq, q\bar{q}
3173
3174
3175 if(PID == 92 || PID == 93 || PID == 94) return kFALSE;
3176 else if(PID == 13 || PID == 28 || PID == 68) return kFALSE;
3177 else if(PID == 11 || PID == 12 || PID == 53) return kTRUE;
3178 else{
3179 AliInfo(Form("unknown Pythia PID %d",PID));
3180 }
3181
3182 return kFALSE;
3183}
3184
3185// ---------------------------------------------------------------------------------
3186
3187Bool_t AliAnalysisTaskJetChem::IsDiffractiveEvent(const Int_t PID){
3188
3189 // Pythia Manual sec. 8.2.1 :
3190 // if Pythia PID = 92,93,94 event is diffractive
3191 // if Pythia PID = 13, 28, 68 hard scattering products are gg, qg
3192 // if Pythia PID = 11, 12, 53 hard scattering products are qq, q\bar{q}
3193
3194 if(PID == -1) return kFALSE;
3195
3196 if(PID == 13 || PID == 28 || PID == 68) return kFALSE;
3197 else if(PID == 11 || PID == 12 || PID == 53) return kFALSE;
3198 else if(PID == 92 || PID == 93 || PID == 94) return kTRUE;
3199 else{
3200 AliInfo(Form("unknown Pythia PID %d",PID));
3201 }
3202
3203 return kFALSE;
3204
3205}
3206
3207
3208// ----------------------------------------------------------------------------------
3209
3210Int_t AliAnalysisTaskJetChem::GetPythiaProcessID(){
3211
3212 // Pythia PID for this event
3213
3214 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
3215 if (!mcHandler) {
3216 //Printf("ERROR: Could not retrieve MC event handler");
3217 return -1;
3218 }
3219
3220 AliMCEvent* mcEvent = mcHandler->MCEvent();
3221 if (!mcEvent) {
3222 Printf("ERROR: Could not retrieve MC event");
3223 return -1;
3224 }
3225
3226 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
3227
3228 if (!mcEvent) {
3229 Printf("ERROR: Could not retrieve pythiaEventHeader");
3230 return -1;
3231 }
3232
3233 Int_t pid = pythiaGenHeader->ProcessType();
3234
3235 return pid;
3236
3237}