JETAN/DEV/AliAnalysisTaskJetCluster.cxx

   1 // **************************************
   2 // Task used for the correction of determiantion of reconstructed jet spectra
   3 // Compares input (gen) and output (rec) jets
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TROOT.h>
  24 #include <TRandom3.h>
  25 #include <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TRefArray.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TF1.h>
  36 #include <TList.h>
  37 #include <TLorentzVector.h>
  38 #include <TClonesArray.h>
  39 #include  "TDatabasePDG.h"
  40
  41 #include "AliAnalysisTaskJetCluster.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliJetFinder.h"
  44 #include "AliJetHeader.h"
  45 #include "AliJetReader.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODExtension.h"
  50 #include "AliAODTrack.h"
  51 #include "AliAODJet.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliJetKineReaderHeader.h"
  58 #include "AliGenCocktailEventHeader.h"
  59 #include "AliInputEventHandler.h"
  60 #include "AliAODJetEventBackground.h"
  61
  62 #include "fastjet/PseudoJet.hh"
  63 #include "fastjet/ClusterSequenceArea.hh"
  64 #include "fastjet/AreaDefinition.hh"
  65 #include "fastjet/JetDefinition.hh"
  66 // get info on how fastjet was configured
  67 #include "fastjet/config.h"
  68
  69 using std::vector;
  70
  71 ClassImp(AliAnalysisTaskJetCluster)
  72
  73 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
  74   //
  75   // Destructor
  76   //
  77
  78   delete fRef;
  79   delete fRandom;
  80
  81   if(fTCAJetsOut)fTCAJetsOut->Delete();
  82   delete fTCAJetsOut;
  83
  84   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
  85   delete fTCAJetsOutRan;
  86
  87   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
  88   delete fTCARandomConesOut;
  89
  90   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
  91   delete fTCARandomConesOutRan;
  92
  93   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
  94   delete fAODJetBackgroundOut;
  95 }
  96
  97 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
  98   AliAnalysisTaskSE(),
  99   fAOD(0x0),
 100   fAODExtension(0x0),
 101   fRef(new TRefArray),
 102   fUseAODTrackInput(kFALSE),
 103   fUseAODMCInput(kFALSE),
 104   fUseBackgroundCalc(kFALSE),
 105   fEventSelection(kFALSE),
 106   fFilterMask(0),
 107   fFilterMaskBestPt(0),
 108   fFilterType(0),
 109   fJetTypes(kJet),
 110   fTrackTypeRec(kTrackUndef),
 111   fTrackTypeGen(kTrackUndef),
 112   fNSkipLeadingRan(0),
 113   fNSkipLeadingCone(0),
 114   fNRandomCones(0),
 115   fAvgTrials(1),
 116   fExternalWeight(1),
 117   fTrackEtaWindow(0.9),
 118   fRecEtaWindow(0.5),
 119   fTrackPtCut(0.),
 120   fJetOutputMinPt(0.150),
 121   fMaxTrackPtInJet(100.),
 122   fJetTriggerPtCut(0),
 123   fVtxZCut(8),
 124   fVtxR2Cut(1),
 125   fCentCutUp(0),
 126   fCentCutLo(0),
 127   fStoreRhoLeadingTrackCorr(kFALSE),
 128   fNonStdBranch(""),
 129   fBackgroundBranch(""),
 130   fNonStdFile(""),
 131   fMomResH1(0x0),
 132   fMomResH2(0x0),
 133   fMomResH3(0x0),
 134   fMomResH1Fit(0x0),
 135   fMomResH2Fit(0x0),
 136   fMomResH3Fit(0x0),
 137   fhEffH1(0x0),
 138   fhEffH2(0x0),
 139   fhEffH3(0x0),
 140   fUseTrPtResolutionSmearing(kFALSE),
 141   fUseDiceEfficiency(kFALSE),
 142   fDiceEfficiencyMinPt(0.),
 143   fUseTrPtResolutionFromOADB(kFALSE),
 144   fUseTrEfficiencyFromOADB(kFALSE),
 145   fPathTrPtResolution(""),
 146   fPathTrEfficiency(""),
 147   fChangeEfficiencyFraction(0.),
 148   fRparam(1.0),
 149   fAlgorithm(fastjet::kt_algorithm),
 150   fStrategy(fastjet::Best),
 151   fRecombScheme(fastjet::BIpt_scheme),
 152   fAreaType(fastjet::active_area),
 153   fGhostArea(0.01),
 154   fActiveAreaRepeats(1),
 155   fGhostEtamax(1.5),
 156   fTCAJetsOut(0x0),
 157   fTCAJetsOutRan(0x0),
 158   fTCARandomConesOut(0x0),
 159   fTCARandomConesOutRan(0x0),
 160   fAODJetBackgroundOut(0x0),
 161   fRandom(0),
 162   fh1Xsec(0x0),
 163   fh1Trials(0x0),
 164   fh1PtHard(0x0),
 165   fh1PtHardNoW(0x0),
 166   fh1PtHardTrials(0x0),
 167   fh1NJetsRec(0x0),
 168   fh1NConstRec(0x0),
 169   fh1NConstLeadingRec(0x0),
 170   fh1PtJetsRecIn(0x0),
 171   fh1PtJetsLeadingRecIn(0x0),
 172   fh1PtJetConstRec(0x0),
 173   fh1PtJetConstLeadingRec(0x0),
 174   fh1PtTracksRecIn(0x0),
 175   fh1PtTracksLeadingRecIn(0x0),
 176   fh1NJetsRecRan(0x0),
 177   fh1NConstRecRan(0x0),
 178   fh1PtJetsLeadingRecInRan(0x0),
 179   fh1NConstLeadingRecRan(0x0),
 180   fh1PtJetsRecInRan(0x0),
 181   fh1PtTracksGenIn(0x0),
 182   fh1Nch(0x0),
 183   fh1CentralityPhySel(0x0),
 184   fh1Centrality(0x0),
 185   fh1CentralitySelect(0x0),
 186   fh1ZPhySel(0x0),
 187   fh1Z(0x0),
 188   fh1ZSelect(0x0),
 189   fh2NRecJetsPt(0x0),
 190   fh2NRecTracksPt(0x0),
 191   fh2NConstPt(0x0),
 192   fh2NConstLeadingPt(0x0),
 193   fh2JetPhiEta(0x0),
 194   fh2LeadingJetPhiEta(0x0),
 195   fh2JetEtaPt(0x0),
 196   fh2LeadingJetEtaPt(0x0),
 197   fh2TrackEtaPt(0x0),
 198   fh2LeadingTrackEtaPt(0x0),
 199   fh2JetsLeadingPhiEta(0x0),
 200   fh2JetsLeadingPhiPt(0x0),
 201   fh2TracksLeadingPhiEta(0x0),
 202   fh2TracksLeadingPhiPt(0x0),
 203   fh2TracksLeadingJetPhiPt(0x0),
 204   fh2JetsLeadingPhiPtW(0x0),
 205   fh2TracksLeadingPhiPtW(0x0),
 206   fh2TracksLeadingJetPhiPtW(0x0),
 207   fh2NRecJetsPtRan(0x0),
 208   fh2NConstPtRan(0x0),
 209   fh2NConstLeadingPtRan(0x0),
 210   fh2PtNch(0x0),
 211   fh2PtNchRan(0x0),
 212   fh2PtNchN(0x0),
 213   fh2PtNchNRan(0x0),
 214   fh2TracksLeadingJetPhiPtRan(0x0),
 215   fh2TracksLeadingJetPhiPtWRan(0x0),
 216   fh2CentvsRho(0x0),
 217   fh2CentvsSigma(0x0),
 218   fh2MultvsRho(0x0),
 219   fh2MultvsSigma(0x0),
 220   fh3CentvsRhoLeadingTrackPtQ1(0x0),
 221   fh3CentvsRhoLeadingTrackPtQ2(0x0),
 222   fh3CentvsRhoLeadingTrackPtQ3(0x0),
 223   fh3CentvsRhoLeadingTrackPtQ4(0x0),
 224   fh3CentvsSigmaLeadingTrackPtQ1(0x0),
 225   fh3CentvsSigmaLeadingTrackPtQ2(0x0),
 226   fh3CentvsSigmaLeadingTrackPtQ3(0x0),
 227   fh3CentvsSigmaLeadingTrackPtQ4(0x0),
 228   fh3MultvsRhoLeadingTrackPtQ1(0x0),
 229   fh3MultvsRhoLeadingTrackPtQ2(0x0),
 230   fh3MultvsRhoLeadingTrackPtQ3(0x0),
 231   fh3MultvsRhoLeadingTrackPtQ4(0x0),
 232   fh3MultvsSigmaLeadingTrackPtQ1(0x0),
 233   fh3MultvsSigmaLeadingTrackPtQ2(0x0),
 234   fh3MultvsSigmaLeadingTrackPtQ3(0x0),
 235   fh3MultvsSigmaLeadingTrackPtQ4(0x0),
 236   fh2PtGenPtSmeared(0x0),
 237   fp1Efficiency(0x0),
 238   fp1PtResolution(0x0),
 239   fHistList(0x0)
 240 {
 241   //
 242   // Constructor
 243   //
 244
 245   for(int i = 0;i<3;i++){
 246     fh1BiARandomCones[i] = 0;
 247     fh1BiARandomConesRan[i] = 0;
 248   }
 249   for(int i = 0;i<kMaxCent;i++){
 250     fh2JetsLeadingPhiPtC[i] = 0;
 251     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 252     fh2TracksLeadingJetPhiPtC[i] = 0;
 253     fh2TracksLeadingJetPhiPtWC[i] = 0;
 254   }
 255 }
 256
 257 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
 258   AliAnalysisTaskSE(name),
 259   fAOD(0x0),
 260   fAODExtension(0x0),
 261   fRef(new TRefArray),
 262   fUseAODTrackInput(kFALSE),
 263   fUseAODMCInput(kFALSE),
 264   fUseBackgroundCalc(kFALSE),
 265   fEventSelection(kFALSE),                                                        fFilterMask(0),
 266   fFilterMaskBestPt(0),
 267   fFilterType(0),
 268   fJetTypes(kJet),
 269   fTrackTypeRec(kTrackUndef),
 270   fTrackTypeGen(kTrackUndef),
 271   fNSkipLeadingRan(0),
 272   fNSkipLeadingCone(0),
 273   fNRandomCones(0),
 274   fAvgTrials(1),
 275   fExternalWeight(1),
 276   fTrackEtaWindow(0.9),
 277   fRecEtaWindow(0.5),
 278   fTrackPtCut(0.),
 279   fJetOutputMinPt(0.150),
 280   fMaxTrackPtInJet(100.),
 281   fJetTriggerPtCut(0),
 282   fVtxZCut(8),
 283   fVtxR2Cut(1),
 284   fCentCutUp(0),
 285   fCentCutLo(0),
 286   fStoreRhoLeadingTrackCorr(kFALSE),
 287   fNonStdBranch(""),
 288   fBackgroundBranch(""),
 289   fNonStdFile(""),
 290   fMomResH1(0x0),
 291   fMomResH2(0x0),
 292   fMomResH3(0x0),
 293   fMomResH1Fit(0x0),
 294   fMomResH2Fit(0x0),
 295   fMomResH3Fit(0x0),
 296   fhEffH1(0x0),
 297   fhEffH2(0x0),
 298   fhEffH3(0x0),
 299   fUseTrPtResolutionSmearing(kFALSE),
 300   fUseDiceEfficiency(kFALSE),
 301   fDiceEfficiencyMinPt(0.),
 302   fUseTrPtResolutionFromOADB(kFALSE),
 303   fUseTrEfficiencyFromOADB(kFALSE),
 304   fPathTrPtResolution(""),
 305   fPathTrEfficiency(""),
 306   fChangeEfficiencyFraction(0.),
 307   fRparam(1.0),
 308   fAlgorithm(fastjet::kt_algorithm),
 309   fStrategy(fastjet::Best),
 310   fRecombScheme(fastjet::BIpt_scheme),
 311   fAreaType(fastjet::active_area),
 312   fGhostArea(0.01),
 313   fActiveAreaRepeats(1),
 314   fGhostEtamax(1.5),
 315   fTCAJetsOut(0x0),
 316   fTCAJetsOutRan(0x0),
 317   fTCARandomConesOut(0x0),
 318   fTCARandomConesOutRan(0x0),
 319   fAODJetBackgroundOut(0x0),
 320   fRandom(0),
 321   fh1Xsec(0x0),
 322   fh1Trials(0x0),
 323   fh1PtHard(0x0),
 324   fh1PtHardNoW(0x0),
 325   fh1PtHardTrials(0x0),
 326   fh1NJetsRec(0x0),
 327   fh1NConstRec(0x0),
 328   fh1NConstLeadingRec(0x0),
 329   fh1PtJetsRecIn(0x0),
 330   fh1PtJetsLeadingRecIn(0x0),
 331   fh1PtJetConstRec(0x0),
 332   fh1PtJetConstLeadingRec(0x0),
 333   fh1PtTracksRecIn(0x0),
 334   fh1PtTracksLeadingRecIn(0x0),
 335   fh1NJetsRecRan(0x0),
 336   fh1NConstRecRan(0x0),
 337   fh1PtJetsLeadingRecInRan(0x0),
 338   fh1NConstLeadingRecRan(0x0),
 339   fh1PtJetsRecInRan(0x0),
 340   fh1PtTracksGenIn(0x0),
 341   fh1Nch(0x0),
 342   fh1CentralityPhySel(0x0),
 343   fh1Centrality(0x0),
 344   fh1CentralitySelect(0x0),
 345   fh1ZPhySel(0x0),
 346   fh1Z(0x0),
 347   fh1ZSelect(0x0),
 348   fh2NRecJetsPt(0x0),
 349   fh2NRecTracksPt(0x0),
 350   fh2NConstPt(0x0),
 351   fh2NConstLeadingPt(0x0),
 352   fh2JetPhiEta(0x0),
 353   fh2LeadingJetPhiEta(0x0),
 354   fh2JetEtaPt(0x0),
 355   fh2LeadingJetEtaPt(0x0),
 356   fh2TrackEtaPt(0x0),
 357   fh2LeadingTrackEtaPt(0x0),
 358   fh2JetsLeadingPhiEta(0x0),
 359   fh2JetsLeadingPhiPt(0x0),
 360   fh2TracksLeadingPhiEta(0x0),
 361   fh2TracksLeadingPhiPt(0x0),
 362   fh2TracksLeadingJetPhiPt(0x0),
 363   fh2JetsLeadingPhiPtW(0x0),
 364   fh2TracksLeadingPhiPtW(0x0),
 365   fh2TracksLeadingJetPhiPtW(0x0),
 366   fh2NRecJetsPtRan(0x0),
 367   fh2NConstPtRan(0x0),
 368   fh2NConstLeadingPtRan(0x0),
 369   fh2PtNch(0x0),
 370   fh2PtNchRan(0x0),
 371   fh2PtNchN(0x0),
 372   fh2PtNchNRan(0x0),
 373   fh2TracksLeadingJetPhiPtRan(0x0),
 374   fh2TracksLeadingJetPhiPtWRan(0x0),
 375   fh2CentvsRho(0x0),
 376   fh2CentvsSigma(0x0),
 377   fh2MultvsRho(0x0),
 378   fh2MultvsSigma(0x0),
 379   fh3CentvsRhoLeadingTrackPtQ1(0x0),
 380   fh3CentvsRhoLeadingTrackPtQ2(0x0),
 381   fh3CentvsRhoLeadingTrackPtQ3(0x0),
 382   fh3CentvsRhoLeadingTrackPtQ4(0x0),
 383   fh3CentvsSigmaLeadingTrackPtQ1(0x0),
 384   fh3CentvsSigmaLeadingTrackPtQ2(0x0),
 385   fh3CentvsSigmaLeadingTrackPtQ3(0x0),
 386   fh3CentvsSigmaLeadingTrackPtQ4(0x0),
 387   fh3MultvsRhoLeadingTrackPtQ1(0x0),
 388   fh3MultvsRhoLeadingTrackPtQ2(0x0),
 389   fh3MultvsRhoLeadingTrackPtQ3(0x0),
 390   fh3MultvsRhoLeadingTrackPtQ4(0x0),
 391   fh3MultvsSigmaLeadingTrackPtQ1(0x0),
 392   fh3MultvsSigmaLeadingTrackPtQ2(0x0),
 393   fh3MultvsSigmaLeadingTrackPtQ3(0x0),
 394   fh3MultvsSigmaLeadingTrackPtQ4(0x0),
 395   fh2PtGenPtSmeared(0x0),
 396   fp1Efficiency(0x0),
 397   fp1PtResolution(0x0),
 398   fHistList(0x0)
 399 {
 400   //
 401   // named ctor
 402   //
 403
 404   for(int i = 0;i<3;i++){
 405     fh1BiARandomCones[i] = 0;
 406     fh1BiARandomConesRan[i] = 0;
 407   }
 408   for(int i = 0;i<kMaxCent;i++){
 409     fh2JetsLeadingPhiPtC[i] = 0;
 410     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 411     fh2TracksLeadingJetPhiPtC[i] = 0;
 412     fh2TracksLeadingJetPhiPtWC[i] = 0;
 413   }
 414  DefineOutput(1, TList::Class());
 415 }
 416
 417
 418
 419 Bool_t AliAnalysisTaskJetCluster::Notify()
 420 {
 421   //
 422   // Implemented Notify() to read the cross sections
 423   // and number of trials from pyxsec.root
 424   //
 425   return kTRUE;
 426 }
 427
 428 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 429 {
 430
 431   //
 432   // Create the output container
 433   //
 434
 435   fRandom = new TRandom3(0);
 436
 437
 438   // Connect the AOD
 439
 440
 441   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 442
 443
 444
 445   if(fNonStdBranch.Length()!=0)
 446     {
 447       // only create the output branch if we have a name
 448       // Create a new branch for jets...
 449       //  -> cleared in the UserExec....
 450       // here we can also have the case that the brnaches are written to a separate file
 451
 452       if(fJetTypes&kJet){
 453         fTCAJetsOut = new TClonesArray("AliAODJet", 0);
 454         fTCAJetsOut->SetName(fNonStdBranch.Data());
 455         AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
 456       }
 457
 458       if(fJetTypes&kJetRan){
 459         fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
 460         fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
 461         if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
 462           fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
 463         AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
 464       }
 465
 466       if(fUseBackgroundCalc){
 467         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
 468           fAODJetBackgroundOut = new AliAODJetEventBackground();
 469           fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
 470           if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
 471             fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
 472
 473           AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
 474         }
 475       }
 476       // create the branch for the random cones with the same R
 477       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
 478       if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
 479         cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
 480
 481       if(fNRandomCones>0){
 482         if(fJetTypes&kRC){
 483           if(!AODEvent()->FindListObject(cName.Data())){
 484             fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
 485             fTCARandomConesOut->SetName(cName.Data());
 486             AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
 487           }
 488         }
 489         // create the branch with the random for the random cones on the random event
 490         if(fJetTypes&kRCRan){
 491           cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 492           if(!AODEvent()->FindListObject(cName.Data())){
 493             fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
 494             fTCARandomConesOutRan->SetName(cName.Data());
 495             AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
 496           }
 497         }
 498       }
 499
 500       if(fNonStdFile.Length()!=0){
 501         //
 502         // case that we have an AOD extension we need to fetch the jets from the extended output
 503         // we identify the extension aod event by looking for the branchname
 504         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 505         // case that we have an AOD extension we need can fetch the background maybe from the extended output
 506         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 507       }
 508     }
 509
 510
 511   if(!fHistList)fHistList = new TList();
 512   fHistList->SetOwner();
 513   PostData(1, fHistList); // post data in any case once
 514
 515   Bool_t oldStatus = TH1::AddDirectoryStatus();
 516   TH1::AddDirectory(kFALSE);
 517
 518   //
 519   //  Histogram
 520
 521   const Int_t nBinPt = 100;
 522   Double_t binLimitsPt[nBinPt+1];
 523   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 524     if(iPt == 0){
 525       binLimitsPt[iPt] = 0.0;
 526     }
 527     else {// 1.0
 528       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
 529     }
 530   }
 531
 532   const Int_t nBinPhi = 90;
 533   Double_t binLimitsPhi[nBinPhi+1];
 534   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 535     if(iPhi==0){
 536       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 537     }
 538     else{
 539       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 540     }
 541   }
 542
 543
 544
 545   const Int_t nBinEta = 40;
 546   Double_t binLimitsEta[nBinEta+1];
 547   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 548     if(iEta==0){
 549       binLimitsEta[iEta] = -2.0;
 550     }
 551     else{
 552       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 553     }
 554   }
 555
 556   const int nChMax = 5000;
 557
 558   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 559   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 560
 561   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 562   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 563
 564
 565   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 566   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 567
 568   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 569   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 570   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 571   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 572
 573
 574   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 575   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 576   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 577
 578   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 579   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 580   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 581   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 582   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 583   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 584   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 585   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 586   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 587   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 588
 589   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 590   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
 591   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 592
 593   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 594   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
 595   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 596
 597   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 598   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 599   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 600   //
 601
 602
 603   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 604   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 605   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 606   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 607
 608   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 609   fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 610   fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 611   fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 612
 613
 614
 615   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 616                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 617   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 618                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 619
 620   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 621                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 622   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 623                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 624
 625   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 626                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 627   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 628                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 629
 630
 631
 632   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 633                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 634   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 635                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 636   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 637                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 638   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 639                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 640   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 641                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 642   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 643                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 644
 645   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 646                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 647   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 648                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 649
 650   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 651                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 652   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 653                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 654
 655   if(fStoreRhoLeadingTrackCorr) {
 656     fh2CentvsRho = new TH2F("fh2CentvsRho","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.);
 657     fh2CentvsSigma = new TH2F("fh2CentvsSigma","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.);
 658     fh2MultvsRho = new TH2F("fh2MultvsRho","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.);
 659     fh2MultvsSigma = new TH2F("fh2MultvsSigma","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.);
 660
 661
 662     fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
 663     fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
 664     fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
 665     fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
 666
 667     fh3CentvsSigmaLeadingTrackPtQ1 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ1","centrality vs background density Q1; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
 668     fh3CentvsSigmaLeadingTrackPtQ2 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ2","centrality vs background density Q2; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
 669     fh3CentvsSigmaLeadingTrackPtQ3 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ3","centrality vs background density Q3; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
 670     fh3CentvsSigmaLeadingTrackPtQ4 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ4","centrality vs background density Q4; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
 671
 672     fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
 673     fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
 674     fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
 675     fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
 676
 677     fh3MultvsSigmaLeadingTrackPtQ1 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
 678     fh3MultvsSigmaLeadingTrackPtQ2 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
 679     fh3MultvsSigmaLeadingTrackPtQ3 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
 680     fh3MultvsSigmaLeadingTrackPtQ4 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
 681
 682     fHistList->Add(fh2CentvsRho);
 683     fHistList->Add(fh2CentvsSigma);
 684     fHistList->Add(fh2MultvsRho);
 685     fHistList->Add(fh2MultvsSigma);
 686
 687     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
 688     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
 689     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
 690     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
 691
 692     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
 693     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
 694     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
 695     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
 696
 697     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
 698     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
 699     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
 700     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
 701
 702     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
 703     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
 704     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
 705     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
 706
 707   }
 708
 709   //Detector level effects histos
 710   fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 711
 712   fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
 713   fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
 714
 715   fHistList->Add(fh2PtGenPtSmeared);
 716   fHistList->Add(fp1Efficiency);
 717   fHistList->Add(fp1PtResolution);
 718
 719   if(fNRandomCones>0&&fUseBackgroundCalc){
 720     for(int i = 0;i<3;i++){
 721       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 722       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 723     }
 724   }
 725
 726   for(int i = 0;i < kMaxCent;i++){
 727     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
 728     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
 729     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
 730     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
 731   }
 732
 733   const Int_t saveLevel = 3; // large save level more histos
 734   if(saveLevel>0){
 735     fHistList->Add(fh1Xsec);
 736     fHistList->Add(fh1Trials);
 737
 738     fHistList->Add(fh1NJetsRec);
 739     fHistList->Add(fh1NConstRec);
 740     fHistList->Add(fh1NConstLeadingRec);
 741     fHistList->Add(fh1PtJetsRecIn);
 742     fHistList->Add(fh1PtJetsLeadingRecIn);
 743     fHistList->Add(fh1PtTracksRecIn);
 744     fHistList->Add(fh1PtTracksLeadingRecIn);
 745     fHistList->Add(fh1PtJetConstRec);
 746     fHistList->Add(fh1PtJetConstLeadingRec);
 747     fHistList->Add(fh1NJetsRecRan);
 748     fHistList->Add(fh1NConstRecRan);
 749     fHistList->Add(fh1PtJetsLeadingRecInRan);
 750     fHistList->Add(fh1NConstLeadingRecRan);
 751     fHistList->Add(fh1PtJetsRecInRan);
 752     fHistList->Add(fh1Nch);
 753     fHistList->Add(fh1Centrality);
 754     fHistList->Add(fh1CentralitySelect);
 755     fHistList->Add(fh1CentralityPhySel);
 756     fHistList->Add(fh1Z);
 757     fHistList->Add(fh1ZSelect);
 758     fHistList->Add(fh1ZPhySel);
 759     if(fNRandomCones>0&&fUseBackgroundCalc){
 760       for(int i = 0;i<3;i++){
 761         fHistList->Add(fh1BiARandomCones[i]);
 762         fHistList->Add(fh1BiARandomConesRan[i]);
 763       }
 764     }
 765   for(int i = 0;i < kMaxCent;i++){
 766     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
 767     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
 768     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
 769     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
 770   }
 771
 772     fHistList->Add(fh2NRecJetsPt);
 773     fHistList->Add(fh2NRecTracksPt);
 774     fHistList->Add(fh2NConstPt);
 775     fHistList->Add(fh2NConstLeadingPt);
 776     fHistList->Add(fh2PtNch);
 777     fHistList->Add(fh2PtNchRan);
 778     fHistList->Add(fh2PtNchN);
 779     fHistList->Add(fh2PtNchNRan);
 780     fHistList->Add(fh2JetPhiEta);
 781     fHistList->Add(fh2LeadingJetPhiEta);
 782     fHistList->Add(fh2JetEtaPt);
 783     fHistList->Add(fh2LeadingJetEtaPt);
 784     fHistList->Add(fh2TrackEtaPt);
 785     fHistList->Add(fh2LeadingTrackEtaPt);
 786     fHistList->Add(fh2JetsLeadingPhiEta );
 787     fHistList->Add(fh2JetsLeadingPhiPt);
 788     fHistList->Add(fh2TracksLeadingPhiEta);
 789     fHistList->Add(fh2TracksLeadingPhiPt);
 790     fHistList->Add(fh2TracksLeadingJetPhiPt);
 791     fHistList->Add(fh2JetsLeadingPhiPtW);
 792     fHistList->Add(fh2TracksLeadingPhiPtW);
 793     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 794     fHistList->Add(fh2NRecJetsPtRan);
 795     fHistList->Add(fh2NConstPtRan);
 796     fHistList->Add(fh2NConstLeadingPtRan);
 797     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 798     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 799   }
 800
 801   // =========== Switch on Sumw2 for all histos ===========
 802   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 803     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 804     if (h1){
 805       h1->Sumw2();
 806       continue;
 807     }
 808     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 809     if(hn)hn->Sumw2();
 810   }
 811   TH1::AddDirectory(oldStatus);
 812 }
 813
 814 void AliAnalysisTaskJetCluster::LocalInit()
 815 {
 816   //
 817   // Initialization
 818   //
 819
 820   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 821
 822   if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
 823   if(fUseTrEfficiencyFromOADB)   LoadTrEfficiencyRootFileFromOADB();
 824
 825
 826   FitMomentumResolution();
 827
 828 }
 829
 830 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 831 {
 832
 833   // handle and reset the output jet branch
 834
 835   if(fTCAJetsOut)fTCAJetsOut->Delete();
 836   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
 837   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
 838   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
 839   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
 840
 841   AliAODJetEventBackground* externalBackground = 0;
 842   if(!externalBackground&&fBackgroundBranch.Length()){
 843     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
 844     if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 845     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 846   }
 847   //
 848   // Execute analysis for current event
 849   //
 850   AliESDEvent *fESD = 0;
 851   if(fUseAODTrackInput){
 852     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 853     if(!fAOD){
 854       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 855       return;
 856     }
 857     // fetch the header
 858   }
 859   else{
 860     //  assume that the AOD is in the general output...
 861     fAOD  = AODEvent();
 862     if(!fAOD){
 863       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 864       return;
 865     }
 866     if(fDebug>0){
 867       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 868     }
 869   }
 870
 871   //Check if information is provided detector level effects
 872   if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
 873   if(!fhEffH1 || !fhEffH2 || !fhEffH3)       fUseDiceEfficiency = kFALSE;
 874
 875   Bool_t selectEvent =  false;
 876   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 877
 878   Float_t cent = 0;
 879   Float_t zVtx  = 0;
 880   Int_t cenClass = -1;
 881   if(fAOD){
 882     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 883     TString vtxTitle(vtxAOD->GetTitle());
 884     zVtx = vtxAOD->GetZ();
 885
 886     cent = fAOD->GetHeader()->GetCentrality();
 887     if(cent<10)cenClass = 0;
 888     else if(cent<30)cenClass = 1;
 889     else if(cent<50)cenClass = 2;
 890     else if(cent<80)cenClass = 3;
 891     if(physicsSelection){
 892       fh1CentralityPhySel->Fill(cent);
 893       fh1ZPhySel->Fill(zVtx);
 894     }
 895
 896     if(fEventSelection){
 897       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 898         Float_t yvtx = vtxAOD->GetY();
 899         Float_t xvtx = vtxAOD->GetX();
 900         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 901         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
 902           if(physicsSelection){
 903             selectEvent = true;
 904           }
 905         }
 906       }
 907       if(fCentCutUp>0){
 908         if(cent<fCentCutLo||cent>fCentCutUp){
 909           selectEvent = false;
 910         }
 911       }
 912     }else{
 913       selectEvent = true;
 914     }
 915   }
 916
 917
 918   if(!selectEvent){
 919     PostData(1, fHistList);
 920     return;
 921   }
 922   fh1Centrality->Fill(cent);
 923   fh1Z->Fill(zVtx);
 924   fh1Trials->Fill("#sum{ntrials}",1);
 925
 926
 927   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 928
 929   // ==== General variables needed
 930
 931
 932
 933   // we simply fetch the tracks/mc particles as a list of AliVParticles
 934
 935   TList recParticles;
 936   TList genParticles;
 937
 938   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 939   Float_t nCh = recParticles.GetEntries();
 940   fh1Nch->Fill(nCh);
 941   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 942   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 943   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 944
 945   // find the jets....
 946
 947   vector<fastjet::PseudoJet> inputParticlesRec;
 948   vector<fastjet::PseudoJet> inputParticlesRecRan;
 949
 950   // Generate the random cones
 951
 952   AliAODJet vTmpRan(1,0,0,1);
 953   for(int i = 0; i < recParticles.GetEntries(); i++){
 954     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 955
 956     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 957     // we take total momentum here
 958
 959     if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
 960       //Add particles to fastjet in case we are not running toy model
 961       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 962       jInp.set_user_index(i);
 963       inputParticlesRec.push_back(jInp);
 964     }
 965     else if(fUseDiceEfficiency) {
 966
 967       // Dice to decide if particle is kept or not - toy  model for efficiency
 968       //
 969       Double_t rnd = fRandom->Uniform(1.);
 970       Double_t pT = vp->Pt();
 971       Double_t eff[3] = {0.};
 972       Double_t pTtmp = pT;
 973       if(pT>10.) pTtmp = 10.;
 974       Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
 975       Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
 976       Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
 977       Int_t cat[3] = {0};
 978       //Sort efficiencies from large to small
 979       if(eff1>eff2 && eff1>eff3) {
 980         eff[0] = eff1;
 981         cat[0] = 1;
 982         if(eff2>eff3) {
 983           eff[1] = eff2;
 984           eff[2] = eff3;
 985           cat[1] = 2;
 986           cat[2] = 3;
 987         } else {
 988           eff[1] = eff3;
 989           eff[2] = eff2;
 990           cat[1] = 3;
 991           cat[2] = 2;
 992         }
 993       }
 994       else if(eff2>eff1 && eff2>eff3) {
 995         eff[0] = eff2;
 996         cat[0] = 2;
 997         if(eff1>eff3) {
 998           eff[1] = eff1;
 999           eff[2] = eff3;
1000           cat[1] = 1;
1001           cat[2] = 3;
1002         } else {
1003           eff[1] = eff3;
1004           eff[2] = eff1;
1005           cat[1] = 3;
1006           cat[2] = 1;
1007         }
1008       }
1009       else if(eff3>eff1 && eff3>eff2) {
1010         eff[0] = eff3;
1011         cat[0] = 3;
1012         if(eff1>eff2) {
1013           eff[1] = eff1;
1014           eff[2] = eff2;
1015           cat[1] = 1;
1016           cat[2] = 2;
1017         } else {
1018           eff[1] = eff2;
1019           eff[2] = eff1;
1020           cat[1] = 2;
1021           cat[2] = 1;
1022         }
1023       }
1024
1025       Double_t sumEff = eff[0]+eff[1]+eff[2];
1026       fp1Efficiency->Fill(vp->Pt(),sumEff);
1027       if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
1028
1029       if(fUseTrPtResolutionSmearing) {
1030         //Smear momentum of generated particle
1031         Double_t smear = 1.;
1032         //Select hybrid track category
1033         if(rnd<=eff[2])
1034           smear = GetMomentumSmearing(cat[2],pT);
1035         else if(rnd<=(eff[2]+eff[1]))
1036           smear = GetMomentumSmearing(cat[1],pT);
1037         else
1038           smear = GetMomentumSmearing(cat[0],pT);
1039
1040         fp1PtResolution->Fill(vp->Pt(),smear);
1041
1042         Double_t sigma = vp->Pt()*smear;
1043         Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
1044
1045         Double_t phi   = vp->Phi();
1046         Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
1047         Double_t pX    = pTrec * TMath::Cos(phi);
1048         Double_t pY    = pTrec * TMath::Sin(phi);
1049         Double_t pZ    = pTrec/TMath::Tan(theta);
1050         Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
1051
1052         fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
1053
1054         fastjet::PseudoJet jInp(pX,pY,pZ,p);
1055         jInp.set_user_index(i);
1056         inputParticlesRec.push_back(jInp);
1057
1058       }
1059       else {
1060         fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1061         jInp.set_user_index(i);
1062         inputParticlesRec.push_back(jInp);
1063
1064       }
1065
1066     }
1067
1068     // the randomized input changes eta and phi, but keeps the p_T
1069     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1070       Double_t pT = vp->Pt();
1071       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
1072       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
1073
1074       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
1075       Double_t pZ = pT/TMath::Tan(theta);
1076
1077       Double_t pX = pT * TMath::Cos(phi);
1078       Double_t pY = pT * TMath::Sin(phi);
1079       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
1080       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1081
1082       jInpRan.set_user_index(i);
1083       inputParticlesRecRan.push_back(jInpRan);
1084       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1085     }
1086
1087     // fill the tref array, only needed when we write out jets
1088     if(fTCAJetsOut){
1089       if(i == 0){
1090         fRef->Delete(); // make sure to delete before placement new...
1091         if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1092           new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1093         }
1094       }
1095       if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp);  //TRefArray does not work with toy model ...
1096     }
1097   }// recparticles
1098
1099   if(inputParticlesRec.size()==0){
1100     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1101     PostData(1, fHistList);
1102     return;
1103   }
1104
1105   // run fast jet
1106   // employ setters for these...
1107
1108
1109   // now create the object that holds info about ghosts
1110   /*
1111   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1112     // reduce CPU time...
1113     fGhostArea = 0.5;
1114     fActiveAreaRepeats = 0;
1115   }
1116   */
1117
1118   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1119   fastjet::AreaType areaType =   fastjet::active_area;
1120   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1121   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1122   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1123
1124   //range where to compute background
1125   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1126   phiMin = 0;
1127   phiMax = 2*TMath::Pi();
1128   rapMax = fGhostEtamax - fRparam;
1129   rapMin = - fGhostEtamax + fRparam;
1130   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1131
1132
1133   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1134   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1135
1136
1137   fh1NJetsRec->Fill(sortedJets.size());
1138
1139  // loop over all jets an fill information, first on is the leading jet
1140
1141   Int_t nRecOver = inclusiveJets.size();
1142   Int_t nRec     = inclusiveJets.size();
1143   if(inclusiveJets.size()>0){
1144     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1145     Double_t area = clustSeq.area(sortedJets[0]);
1146     leadingJet.SetEffArea(area,0);
1147     Float_t pt = leadingJet.Pt();
1148     Int_t nAodOutJets = 0;
1149     Int_t nAodOutTracks = 0;
1150     AliAODJet *aodOutJet = 0;
1151
1152     Int_t iCount = 0;
1153     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1154       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1155       while(pt<ptCut&&iCount<nRec){
1156         nRecOver--;
1157         iCount++;
1158         if(iCount<nRec){
1159           pt = sortedJets[iCount].perp();
1160         }
1161       }
1162       if(nRecOver<=0)break;
1163       fh2NRecJetsPt->Fill(ptCut,nRecOver);
1164     }
1165     Float_t phi = leadingJet.Phi();
1166     if(phi<0)phi+=TMath::Pi()*2.;
1167     Float_t eta = leadingJet.Eta();
1168     Float_t pTback = 0;
1169     if(externalBackground){
1170       // carefull has to be filled in a task before
1171       // todo, ReArrange to the botom
1172      pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1173     }
1174     pt = leadingJet.Pt() - pTback;
1175     // correlation of leading jet with tracks
1176     TIterator *recIter = recParticles.MakeIterator();
1177     recIter->Reset();
1178     AliVParticle *tmpRecTrack = 0;
1179     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1180       Float_t tmpPt = tmpRecTrack->Pt();
1181       // correlation
1182       Float_t tmpPhi =  tmpRecTrack->Phi();
1183       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1184       Float_t dPhi = phi - tmpPhi;
1185       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1186       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1187       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1188       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1189       if(cenClass>=0){
1190         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1191         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1192       }
1193
1194     }
1195
1196
1197     TLorentzVector vecareab;
1198     for(int j = 0; j < nRec;j++){
1199       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1200       aodOutJet = 0;
1201       nAodOutTracks = 0;
1202       Float_t tmpPt = tmpRec.Pt();
1203
1204       if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1205         aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1206         aodOutJet->GetRefTracks()->Clear();
1207         Double_t area1 = clustSeq.area(sortedJets[j]);
1208         aodOutJet->SetEffArea(area1,0);
1209         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1210         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1211         aodOutJet->SetVectorAreaCharged(&vecareab);
1212       }
1213
1214
1215       Float_t tmpPtBack = 0;
1216       if(externalBackground){
1217         // carefull has to be filled in a task before
1218        // todo, ReArrange to the botom
1219         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1220       }
1221       tmpPt = tmpPt - tmpPtBack;
1222       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1223
1224       fh1PtJetsRecIn->Fill(tmpPt);
1225       // Fill Spectra with constituentsemacs
1226       const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1227
1228       fh1NConstRec->Fill(constituents.size());
1229       fh2PtNch->Fill(nCh,tmpPt);
1230       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1231       fh2NConstPt->Fill(tmpPt,constituents.size());
1232       // loop over constiutents and fill spectrum
1233
1234       AliVParticle *partLead = 0;
1235       Float_t ptLead = -1;
1236
1237       for(unsigned int ic = 0; ic < constituents.size();ic++){
1238         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1239         if(!part) continue;
1240         fh1PtJetConstRec->Fill(part->Pt());
1241         if(aodOutJet){
1242           if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1243           if(part->Pt()>fMaxTrackPtInJet){
1244             aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1245           }
1246         }
1247         if(part->Pt()>ptLead){
1248           ptLead = part->Pt();
1249           partLead = part;
1250         }
1251         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1252       }
1253
1254       AliAODTrack *aodT = 0;
1255       if(partLead){
1256         if(aodOutJet){
1257           //set pT of leading constituent of jet
1258           aodOutJet->SetPtLeading(partLead->Pt());
1259           aodT = dynamic_cast<AliAODTrack*>(partLead);
1260           if(aodT){
1261             if(aodT->TestFilterBit(fFilterMaskBestPt)){
1262               aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1263             }
1264           }
1265         }
1266       }
1267
1268      // correlation
1269      Float_t tmpPhi =  tmpRec.Phi();
1270      Float_t tmpEta =  tmpRec.Eta();
1271      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1272      if(j==0){
1273        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1274        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1275        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1276        fh1NConstLeadingRec->Fill(constituents.size());
1277        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1278        continue;
1279      }
1280      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1281      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1282      Float_t dPhi = phi - tmpPhi;
1283      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1284      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1285      Float_t dEta = eta - tmpRec.Eta();
1286      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1287      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1288      if(cenClass>=0){
1289        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1290        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1291      }
1292      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1293     }// loop over reconstructed jets
1294    delete recIter;
1295
1296
1297
1298    // Add the random cones...
1299    if(fNRandomCones>0&&fTCARandomConesOut){
1300      // create a random jet within the acceptance
1301      Double_t etaMax = fTrackEtaWindow - fRparam;
1302      Int_t nCone = 0;
1303      Int_t nConeRan = 0;
1304      Double_t pTC = 1; // small number
1305      for(int ir = 0;ir < fNRandomCones;ir++){
1306        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1307        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1308        // massless jet
1309        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1310        Double_t pZC = pTC/TMath::Tan(thetaC);
1311        Double_t pXC = pTC * TMath::Cos(phiC);
1312        Double_t pYC = pTC * TMath::Sin(phiC);
1313        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1314        AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1315        bool skip = false;
1316        for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1317          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1318          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1319            skip = true;
1320            break;
1321          }
1322        }
1323        // test for overlap with previous cones to avoid double counting
1324        for(int iic = 0;iic<ir;iic++){
1325          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1326          if(iicone){
1327            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1328              skip = true;
1329              break;
1330            }
1331          }
1332        }
1333        if(skip)continue;
1334        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1335        tmpRecC.SetPtLeading(-1.);
1336        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1337        if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1338      }// loop over random cones creation
1339
1340
1341      // loop over the reconstructed particles and add up the pT in the random cones
1342      // maybe better to loop over randomized particles not in the real jets...
1343      // but this by definition brings dow average energy in the whole  event
1344      AliAODJet vTmpRanR(1,0,0,1);
1345      for(int i = 0; i < recParticles.GetEntries(); i++){
1346        AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1347        if(fTCARandomConesOut){
1348          for(int ir = 0;ir < fNRandomCones;ir++){
1349            AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1350            if(jC&&jC->DeltaR(vp)<fRparam){
1351              if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1352              jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1353              if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1354            }
1355          }
1356        }// add up energy in cone
1357
1358        // the randomized input changes eta and phi, but keeps the p_T
1359        if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1360          Double_t pTR = vp->Pt();
1361          Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1362          Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1363
1364          Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1365          Double_t pZR = pTR/TMath::Tan(thetaR);
1366
1367          Double_t pXR = pTR * TMath::Cos(phiR);
1368          Double_t pYR = pTR * TMath::Sin(phiR);
1369          Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
1370          vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1371          if(fTCARandomConesOutRan){
1372            for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1373              AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1374              if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1375                if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1376                jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1377                if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1378              }
1379            }
1380          }
1381        }
1382      }// loop over recparticles
1383
1384      Float_t jetArea = fRparam*fRparam*TMath::Pi();
1385      if(fTCARandomConesOut){
1386        for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1387          // rescale the momentum vectors for the random cones
1388
1389          AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1390          if(rC){
1391            Double_t etaC = rC->Eta();
1392            Double_t phiC = rC->Phi();
1393            // massless jet, unit vector
1394            pTC = rC->ChargedBgEnergy();
1395            if(pTC<=0)pTC = 0.001; // for almost empty events
1396            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1397            Double_t pZC = pTC/TMath::Tan(thetaC);
1398            Double_t pXC = pTC * TMath::Cos(phiC);
1399            Double_t pYC = pTC * TMath::Sin(phiC);
1400            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1401            rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1402            rC->SetBgEnergy(0,0);
1403            rC->SetEffArea(jetArea,0);
1404          }
1405        }
1406      }
1407      if(fTCARandomConesOutRan){
1408        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1409          AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1410          // same wit random
1411          if(rC){
1412            Double_t etaC = rC->Eta();
1413            Double_t phiC = rC->Phi();
1414            // massless jet, unit vector
1415            pTC = rC->ChargedBgEnergy();
1416            if(pTC<=0)pTC = 0.001;// for almost empty events
1417            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1418            Double_t pZC = pTC/TMath::Tan(thetaC);
1419            Double_t pXC = pTC * TMath::Cos(phiC);
1420            Double_t pYC = pTC * TMath::Sin(phiC);
1421            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1422            rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1423            rC->SetBgEnergy(0,0);
1424            rC->SetEffArea(jetArea,0);
1425          }
1426        }
1427      }
1428    }// if(fNRandomCones
1429
1430    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1431
1432    if(fAODJetBackgroundOut){
1433      vector<fastjet::PseudoJet> jets2=sortedJets;
1434      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1435
1436      Double_t bkg1=0;
1437      Double_t sigma1=0.;
1438      Double_t meanarea1=0.;
1439      Double_t bkg2=0;
1440      Double_t sigma2=0.;
1441      Double_t meanarea2=0.;
1442
1443      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1444      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1445
1446      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1447      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1448
1449      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1450      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1451      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1452      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1453
1454    }
1455   }
1456
1457   if(fStoreRhoLeadingTrackCorr) {
1458     vector<fastjet::PseudoJet> jets3=sortedJets;
1459     if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2);
1460
1461     Double_t bkg2=0;
1462     Double_t sigma2=0.;
1463     Double_t meanarea2=0.;
1464
1465     clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1466     fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1467
1468     fh2CentvsRho->Fill(cent,bkg2);
1469     fh2CentvsSigma->Fill(cent,sigma2);
1470     fh2MultvsRho->Fill(nCh,bkg2);
1471     fh2MultvsSigma->Fill(nCh,sigma2);
1472
1473
1474     //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1475     //exclude 2 leading jets from event
1476     //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R  (Near side to leading track)
1477     //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1478     //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R)  (Away side to leading track)
1479     //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R)   (Orthogonal to leading track)
1480
1481     //Assuming R=0.2 for background clusters
1482
1483     Double_t bkg2Q[4]      = {0.};
1484     Double_t sigma2Q[4]    = {0.};
1485     Double_t meanarea2Q[4] = {0.};
1486
1487     //Get phi of leading track in event
1488     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1489     Float_t phiLeadingTrack = leading->Phi();
1490     Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1491
1492     //Quadrant1 - near side
1493     phiMin = phiLeadingTrack - phiStep;
1494     phiMax = phiLeadingTrack + phiStep;
1495     fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1496     clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1497
1498     //Quadrant2 - orthogonal
1499     Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1500     if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1501     phiMin = phiQ2 - phiStep;
1502     phiMax = phiQ2 + phiStep;
1503     fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1504     clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1505
1506     //Quadrant3 - away side
1507     Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1508     if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1509     phiMin = phiQ3 - phiStep;
1510     phiMax = phiQ3 + phiStep;
1511     fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1512     clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1513
1514     //Quadrant4 - othogonal
1515     Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1516     if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1517     phiMin = phiQ4 - phiStep;
1518     phiMax = phiQ4 + phiStep;
1519     fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1520     clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1521
1522     fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1523     fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1524     fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1525     fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1526
1527     fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1528     fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1529     fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1530     fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1531
1532     fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1533     fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1534     fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1535     fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1536
1537     fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1538     fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1539     fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1540     fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1541
1542   }
1543
1544
1545
1546
1547
1548   // fill track information
1549   Int_t nTrackOver = recParticles.GetSize();
1550   // do the same for tracks and jets
1551
1552   if(nTrackOver>0){
1553    TIterator *recIter = recParticles.MakeIterator();
1554    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1555    Float_t pt = tmpRec->Pt();
1556
1557    //    Printf("Leading track p_t %3.3E",pt);
1558    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1559      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1560      while(pt<ptCut&&tmpRec){
1561        nTrackOver--;
1562        tmpRec = (AliVParticle*)(recIter->Next());
1563        if(tmpRec){
1564          pt = tmpRec->Pt();
1565        }
1566      }
1567      if(nTrackOver<=0)break;
1568      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1569    }
1570
1571    recIter->Reset();
1572    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1573    Float_t phi = leading->Phi();
1574    if(phi<0)phi+=TMath::Pi()*2.;
1575    Float_t eta = leading->Eta();
1576    pt  = leading->Pt();
1577    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1578      Float_t tmpPt = tmpRec->Pt();
1579      Float_t tmpEta = tmpRec->Eta();
1580      fh1PtTracksRecIn->Fill(tmpPt);
1581      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1582      if(tmpRec==leading){
1583        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1584        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1585        continue;
1586      }
1587       // correlation
1588      Float_t tmpPhi =  tmpRec->Phi();
1589
1590      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1591      Float_t dPhi = phi - tmpPhi;
1592      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1593      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1594      Float_t dEta = eta - tmpRec->Eta();
1595      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1596      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1597      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1598    }
1599    delete recIter;
1600  }
1601
1602  // find the random jets
1603
1604  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1605
1606  // fill the jet information from random track
1607  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1608  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1609
1610   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1611
1612  // loop over all jets an fill information, first on is the leading jet
1613
1614  Int_t nRecOverRan = inclusiveJetsRan.size();
1615  Int_t nRecRan     = inclusiveJetsRan.size();
1616
1617  if(inclusiveJetsRan.size()>0){
1618    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1619    Float_t pt = leadingJet.Pt();
1620
1621    Int_t iCount = 0;
1622    TLorentzVector vecarearanb;
1623
1624    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1625      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1626       while(pt<ptCut&&iCount<nRecRan){
1627         nRecOverRan--;
1628         iCount++;
1629         if(iCount<nRecRan){
1630           pt = sortedJetsRan[iCount].perp();
1631         }
1632       }
1633       if(nRecOverRan<=0)break;
1634       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1635     }
1636     Float_t phi = leadingJet.Phi();
1637     if(phi<0)phi+=TMath::Pi()*2.;
1638     pt = leadingJet.Pt();
1639
1640     // correlation of leading jet with random tracks
1641
1642     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1643       {
1644         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1645         // correlation
1646         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1647         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1648         Float_t dPhi = phi - tmpPhi;
1649         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1650         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1651         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1652         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1653       }
1654
1655     Int_t nAodOutJetsRan = 0;
1656      AliAODJet *aodOutJetRan = 0;
1657     for(int j = 0; j < nRecRan;j++){
1658       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1659       Float_t tmpPt = tmpRec.Pt();
1660       fh1PtJetsRecInRan->Fill(tmpPt);
1661       // Fill Spectra with constituents
1662       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1663       fh1NConstRecRan->Fill(constituents.size());
1664       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1665       fh2PtNchRan->Fill(nCh,tmpPt);
1666       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1667
1668
1669      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1670        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1671        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1672        aodOutJetRan->GetRefTracks()->Clear();
1673        aodOutJetRan->SetEffArea(arearan,0);
1674        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1675        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1676        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1677
1678      }
1679
1680       // correlation
1681       Float_t tmpPhi =  tmpRec.Phi();
1682       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1683
1684       if(j==0){
1685         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1686         fh1NConstLeadingRecRan->Fill(constituents.size());
1687         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1688         continue;
1689       }
1690     }
1691
1692
1693     if(fAODJetBackgroundOut){
1694      Double_t bkg3=0.;
1695      Double_t sigma3=0.;
1696      Double_t meanarea3=0.;
1697      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1698      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1699      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1700      /*
1701      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1702      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1703      */
1704     }
1705
1706
1707
1708  }
1709
1710
1711  // do the event selection if activated
1712  if(fJetTriggerPtCut>0){
1713    bool select = false;
1714    Float_t minPt = fJetTriggerPtCut;
1715    /*
1716    // hard coded for now ...
1717    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1718    if(cent<10)minPt = 50;
1719    else if(cent<30)minPt = 42;
1720    else if(cent<50)minPt = 28;
1721    else if(cent<80)minPt = 18;
1722    */
1723    float rho = 0;
1724    if(externalBackground)rho = externalBackground->GetBackground(2);
1725    if(fTCAJetsOut){
1726      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1727        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1728        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1729        if(ptSub>=minPt){
1730          select = true;
1731          break;
1732        }
1733      }
1734    }
1735
1736    if(select){
1737      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1738      fh1CentralitySelect->Fill(cent);
1739      fh1ZSelect->Fill(zVtx);
1740      aodH->SetFillAOD(kTRUE);
1741    }
1742  }
1743  if (fDebug > 2){
1744    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1745    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1746    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1747    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1748  }
1749  PostData(1, fHistList);
1750 }
1751
1752 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1753 {
1754   //
1755   // Terminate analysis
1756   //
1757     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1758
1759     if(fMomResH1Fit) delete fMomResH1Fit;
1760     if(fMomResH2Fit) delete fMomResH2Fit;
1761     if(fMomResH3Fit) delete fMomResH3Fit;
1762
1763 }
1764
1765
1766 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1767
1768   //
1769   // get list of tracks/particles for different types
1770   //
1771
1772   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1773
1774   Int_t iCount = 0;
1775   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1776     if(type!=kTrackAODextraonly) {
1777       AliAODEvent *aod = 0;
1778       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1779       else aod = AODEvent();
1780       if(!aod){
1781         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1782         return iCount;
1783       }
1784
1785       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1786         AliAODTrack *tr = aod->GetTrack(it);
1787         Bool_t bGood = false;
1788         if(fFilterType == 0)bGood = true;
1789         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1790         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1791         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1792           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1793           continue;
1794         }
1795         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1796           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1797           continue;
1798         }
1799         if(tr->Pt()<fTrackPtCut){
1800           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1801           continue;
1802         }
1803         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1804         list->Add(tr);
1805         iCount++;
1806       }
1807     }
1808     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1809       AliAODEvent *aod = 0;
1810       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1811       else aod = AODEvent();
1812
1813       if(!aod){
1814         return iCount;
1815       }
1816       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1817       if(!aodExtraTracks)return iCount;
1818       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1819         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1820         if (!track) continue;
1821
1822         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1823         if(!trackAOD)continue;
1824         Bool_t bGood = false;
1825         if(fFilterType == 0)bGood = true;
1826         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1827         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1828         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1829         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1830         if(trackAOD->Pt()<fTrackPtCut) continue;
1831         list->Add(trackAOD);
1832         iCount++;
1833       }
1834     }
1835   }
1836   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1837     AliMCEvent* mcEvent = MCEvent();
1838     if(!mcEvent)return iCount;
1839     // we want to have alivpartilces so use get track
1840     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1841       if(!mcEvent->IsPhysicalPrimary(it))continue;
1842       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1843       if(type == kTrackKineAll){
1844         if(part->Pt()<fTrackPtCut)continue;
1845         list->Add(part);
1846         iCount++;
1847       }
1848       else if(type == kTrackKineCharged){
1849         if(part->Particle()->GetPDG()->Charge()==0)continue;
1850         if(part->Pt()<fTrackPtCut)continue;
1851         list->Add(part);
1852         iCount++;
1853       }
1854     }
1855   }
1856   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1857     AliAODEvent *aod = 0;
1858     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1859     else aod = AODEvent();
1860     if(!aod)return iCount;
1861     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1862     if(!tca)return iCount;
1863     for(int it = 0;it < tca->GetEntriesFast();++it){
1864       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1865       if(!part->IsPhysicalPrimary())continue;
1866       if(type == kTrackAODMCAll){
1867         if(part->Pt()<fTrackPtCut)continue;
1868         list->Add(part);
1869         iCount++;
1870       }
1871       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1872         if(part->Charge()==0)continue;
1873         if(part->Pt()<fTrackPtCut)continue;
1874         if(kTrackAODMCCharged){
1875           list->Add(part);
1876         }
1877         else {
1878           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1879           list->Add(part);
1880         }
1881         iCount++;
1882       }
1883     }
1884   }// AODMCparticle
1885   list->Sort();
1886   return iCount;
1887 }
1888
1889 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1890
1891   TFile *f = TFile::Open(fPathTrPtResolution.Data());
1892
1893   if(!f)return;
1894
1895   TProfile *fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1896   TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1897   TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1898
1899   SetSmearResolution(kTRUE);
1900   SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
1901
1902
1903 }
1904
1905 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1906
1907   TFile *f = TFile::Open(fPathTrEfficiency.Data());
1908   if(!f)return;
1909
1910   TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1911   TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1912   TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1913
1914   SetDiceEfficiency(kTRUE);
1915
1916   if(fChangeEfficiencyFraction>0.) {
1917
1918     TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1919
1920     for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1921       Double_t content = hEffPosGlobSt->GetBinContent(bin);
1922       hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1923     }
1924
1925     SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1926
1927   }
1928   else
1929     SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1930
1931 }
1932
1933 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1934
1935   //
1936   // set mom res profiles
1937   //
1938
1939   if(fMomResH1) delete fMomResH1;
1940   if(fMomResH2) delete fMomResH2;
1941   if(fMomResH3) delete fMomResH3;
1942
1943   fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
1944   fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
1945   fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
1946
1947 }
1948
1949 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1950   //
1951   // set tracking efficiency histos
1952   //
1953
1954   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1955   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1956   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1957 }
1958
1959 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1960
1961   //
1962   // Get smearing on generated momentum
1963   //
1964
1965   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1966
1967   TProfile *fMomRes = 0x0;
1968   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1969   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1970   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1971
1972   if(!fMomRes) {
1973     return 0.;
1974   }
1975
1976
1977   Double_t smear = 0.;
1978
1979   if(pt>20.) {
1980     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1981     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1982     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1983   }
1984   else {
1985
1986     Int_t bin = fMomRes->FindBin(pt);
1987
1988     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1989
1990   }
1991
1992   if(fMomRes) delete fMomRes;
1993
1994   return smear;
1995 }
1996
1997 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1998   //
1999   // Fit linear function on momentum resolution at high pT
2000   //
2001
2002   if(!fMomResH1Fit && fMomResH1) {
2003     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2004     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2005     fMomResH1Fit ->SetRange(5.,100.);
2006   }
2007
2008   if(!fMomResH2Fit && fMomResH2) {
2009     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2010     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2011     fMomResH2Fit ->SetRange(5.,100.);
2012   }
2013
2014   if(!fMomResH3Fit && fMomResH3) {
2015     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2016     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2017     fMomResH3Fit ->SetRange(5.,100.);
2018   }
2019
2020 }
2021
2022 /*
2023 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2024   for(int i = 0; i < particles.GetEntries(); i++){
2025     AliVParticle *vp = (AliVParticle*)particles.At(i);
2026     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2027     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2028     jInp.set_user_index(i);
2029     inputParticles.push_back(jInp);
2030   }
2031
2032   return 0;
2033
2034 }
2035 */