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