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