PWG4/JetTasks/AliAnalysisTaskJetChem.cxx

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