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