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