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