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