JETAN/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
1433
1434
1435
1436    if(fAODJetBackgroundOut){
1437      vector<fastjet::PseudoJet> jets2=sortedJets;
1438      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1439      Double_t bkg1=0;
1440      Double_t sigma1=0.;
1441      Double_t meanarea1=0.;
1442      Double_t bkg2=0;
1443      Double_t sigma2=0.;
1444      Double_t meanarea2=0.;
1445
1446      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1447      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1448
1449      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1450      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1451
1452      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1453      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1454      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1455      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1456
1457      if(fStoreRhoLeadingTrackCorr) {
1458        fh2CentvsRho->Fill(cent,bkg2);
1459        fh2CentvsSigma->Fill(cent,sigma2);
1460        fh2MultvsRho->Fill(nCh,bkg2);
1461        fh2MultvsSigma->Fill(nCh,sigma2);
1462
1463
1464        //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1465        //exclude 2 leading jets from event
1466        //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R  (Near side to leading track)
1467        //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1468        //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R)  (Away side to leading track)
1469        //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R)   (Orthogonal to leading track)
1470
1471        //Assuming R=0.2 for background clusters
1472
1473        Double_t bkg2Q[4]      = {0.};
1474        Double_t sigma2Q[4]    = {0.};
1475        Double_t meanarea2Q[4] = {0.};
1476
1477        //Get phi of leading track in event
1478        AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1479        Float_t phiLeadingTrack = leading->Phi();
1480        Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1481
1482        //Quadrant1 - near side
1483        phiMin = phiLeadingTrack - phiStep;
1484        phiMax = phiLeadingTrack + phiStep;
1485        fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1486        clustSeq.get_median_rho_and_sigma(jets2, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1487
1488        //Quadrant2 - orthogonal
1489        Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1490        if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1491        phiMin = phiQ2 - phiStep;
1492        phiMax = phiQ2 + phiStep;
1493        fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1494        clustSeq.get_median_rho_and_sigma(jets2, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1495
1496        //Quadrant3 - away side
1497        Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1498        if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1499        phiMin = phiQ3 - phiStep;
1500        phiMax = phiQ3 + phiStep;
1501        fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1502        clustSeq.get_median_rho_and_sigma(jets2, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1503
1504        //Quadrant4 - othogonal
1505        Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1506        if(phiQ3>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1507        phiMin = phiQ4 - phiStep;
1508        phiMax = phiQ4 + phiStep;
1509        fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1510        clustSeq.get_median_rho_and_sigma(jets2, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1511
1512        fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1513        fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1514        fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1515        fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1516
1517        fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1518        fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1519        fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1520        fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1521
1522        fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1523        fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1524        fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1525        fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1526
1527        fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1528        fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1529        fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1530        fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1531
1532      }
1533
1534
1535    }
1536   }
1537
1538
1539
1540
1541
1542   // fill track information
1543   Int_t nTrackOver = recParticles.GetSize();
1544   // do the same for tracks and jets
1545
1546   if(nTrackOver>0){
1547    TIterator *recIter = recParticles.MakeIterator();
1548    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1549    Float_t pt = tmpRec->Pt();
1550
1551    //    Printf("Leading track p_t %3.3E",pt);
1552    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1553      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1554      while(pt<ptCut&&tmpRec){
1555        nTrackOver--;
1556        tmpRec = (AliVParticle*)(recIter->Next());
1557        if(tmpRec){
1558          pt = tmpRec->Pt();
1559        }
1560      }
1561      if(nTrackOver<=0)break;
1562      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1563    }
1564
1565    recIter->Reset();
1566    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1567    Float_t phi = leading->Phi();
1568    if(phi<0)phi+=TMath::Pi()*2.;
1569    Float_t eta = leading->Eta();
1570    pt  = leading->Pt();
1571    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1572      Float_t tmpPt = tmpRec->Pt();
1573      Float_t tmpEta = tmpRec->Eta();
1574      fh1PtTracksRecIn->Fill(tmpPt);
1575      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1576      if(tmpRec==leading){
1577        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1578        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1579        continue;
1580      }
1581       // correlation
1582      Float_t tmpPhi =  tmpRec->Phi();
1583
1584      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1585      Float_t dPhi = phi - tmpPhi;
1586      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1587      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1588      Float_t dEta = eta - tmpRec->Eta();
1589      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1590      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1591      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1592    }
1593    delete recIter;
1594  }
1595
1596  // find the random jets
1597
1598  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1599
1600  // fill the jet information from random track
1601  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1602  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1603
1604   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1605
1606  // loop over all jets an fill information, first on is the leading jet
1607
1608  Int_t nRecOverRan = inclusiveJetsRan.size();
1609  Int_t nRecRan     = inclusiveJetsRan.size();
1610
1611  if(inclusiveJetsRan.size()>0){
1612    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1613    Float_t pt = leadingJet.Pt();
1614
1615    Int_t iCount = 0;
1616    TLorentzVector vecarearanb;
1617
1618    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1619      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1620       while(pt<ptCut&&iCount<nRecRan){
1621         nRecOverRan--;
1622         iCount++;
1623         if(iCount<nRecRan){
1624           pt = sortedJetsRan[iCount].perp();
1625         }
1626       }
1627       if(nRecOverRan<=0)break;
1628       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1629     }
1630     Float_t phi = leadingJet.Phi();
1631     if(phi<0)phi+=TMath::Pi()*2.;
1632     pt = leadingJet.Pt();
1633
1634     // correlation of leading jet with random tracks
1635
1636     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1637       {
1638         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1639         // correlation
1640         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1641         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1642         Float_t dPhi = phi - tmpPhi;
1643         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1644         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1645         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1646         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1647       }
1648
1649     Int_t nAodOutJetsRan = 0;
1650      AliAODJet *aodOutJetRan = 0;
1651     for(int j = 0; j < nRecRan;j++){
1652       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1653       Float_t tmpPt = tmpRec.Pt();
1654       fh1PtJetsRecInRan->Fill(tmpPt);
1655       // Fill Spectra with constituents
1656       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1657       fh1NConstRecRan->Fill(constituents.size());
1658       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1659       fh2PtNchRan->Fill(nCh,tmpPt);
1660       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1661
1662
1663      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1664        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1665        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1666        aodOutJetRan->GetRefTracks()->Clear();
1667        aodOutJetRan->SetEffArea(arearan,0);
1668        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1669        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1670        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1671
1672      }
1673
1674       // correlation
1675       Float_t tmpPhi =  tmpRec.Phi();
1676       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1677
1678       if(j==0){
1679         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1680         fh1NConstLeadingRecRan->Fill(constituents.size());
1681         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1682         continue;
1683       }
1684     }
1685
1686
1687     if(fAODJetBackgroundOut){
1688      Double_t bkg3=0.;
1689      Double_t sigma3=0.;
1690      Double_t meanarea3=0.;
1691      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1692      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1693      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1694      /*
1695      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1696      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1697      */
1698     }
1699
1700
1701
1702  }
1703
1704
1705  // do the event selection if activated
1706  if(fJetTriggerPtCut>0){
1707    bool select = false;
1708    Float_t minPt = fJetTriggerPtCut;
1709    /*
1710    // hard coded for now ...
1711    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1712    if(cent<10)minPt = 50;
1713    else if(cent<30)minPt = 42;
1714    else if(cent<50)minPt = 28;
1715    else if(cent<80)minPt = 18;
1716    */
1717    float rho = 0;
1718    if(externalBackground)rho = externalBackground->GetBackground(2);
1719    if(fTCAJetsOut){
1720      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1721        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1722        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1723        if(ptSub>=minPt){
1724          select = true;
1725          break;
1726        }
1727      }
1728    }
1729
1730    if(select){
1731      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1732      fh1CentralitySelect->Fill(cent);
1733      fh1ZSelect->Fill(zVtx);
1734      aodH->SetFillAOD(kTRUE);
1735    }
1736  }
1737  if (fDebug > 2){
1738    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1739    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1740    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1741    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1742  }
1743  PostData(1, fHistList);
1744 }
1745
1746 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1747 {
1748   //
1749   // Terminate analysis
1750   //
1751     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1752
1753     if(fMomResH1Fit) delete fMomResH1Fit;
1754     if(fMomResH2Fit) delete fMomResH2Fit;
1755     if(fMomResH3Fit) delete fMomResH3Fit;
1756
1757 }
1758
1759
1760 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1761
1762   //
1763   // get list of tracks/particles for different types
1764   //
1765
1766   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1767
1768   Int_t iCount = 0;
1769   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1770     if(type!=kTrackAODextraonly) {
1771       AliAODEvent *aod = 0;
1772       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1773       else aod = AODEvent();
1774       if(!aod){
1775         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1776         return iCount;
1777       }
1778
1779       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1780         AliAODTrack *tr = aod->GetTrack(it);
1781         Bool_t bGood = false;
1782         if(fFilterType == 0)bGood = true;
1783         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1784         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1785         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1786           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1787           continue;
1788         }
1789         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1790           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1791           continue;
1792         }
1793         if(tr->Pt()<fTrackPtCut){
1794           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1795           continue;
1796         }
1797         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1798         list->Add(tr);
1799         iCount++;
1800       }
1801     }
1802     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1803       AliAODEvent *aod = 0;
1804       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1805       else aod = AODEvent();
1806
1807       if(!aod){
1808         return iCount;
1809       }
1810       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1811       if(!aodExtraTracks)return iCount;
1812       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1813         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1814         if (!track) continue;
1815
1816         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1817         if(!trackAOD)continue;
1818         Bool_t bGood = false;
1819         if(fFilterType == 0)bGood = true;
1820         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1821         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1822         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1823         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1824         if(trackAOD->Pt()<fTrackPtCut) continue;
1825         list->Add(trackAOD);
1826         iCount++;
1827       }
1828     }
1829   }
1830   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1831     AliMCEvent* mcEvent = MCEvent();
1832     if(!mcEvent)return iCount;
1833     // we want to have alivpartilces so use get track
1834     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1835       if(!mcEvent->IsPhysicalPrimary(it))continue;
1836       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1837       if(type == kTrackKineAll){
1838         if(part->Pt()<fTrackPtCut)continue;
1839         list->Add(part);
1840         iCount++;
1841       }
1842       else if(type == kTrackKineCharged){
1843         if(part->Particle()->GetPDG()->Charge()==0)continue;
1844         if(part->Pt()<fTrackPtCut)continue;
1845         list->Add(part);
1846         iCount++;
1847       }
1848     }
1849   }
1850   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1851     AliAODEvent *aod = 0;
1852     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1853     else aod = AODEvent();
1854     if(!aod)return iCount;
1855     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1856     if(!tca)return iCount;
1857     for(int it = 0;it < tca->GetEntriesFast();++it){
1858       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1859       if(!part->IsPhysicalPrimary())continue;
1860       if(type == kTrackAODMCAll){
1861         if(part->Pt()<fTrackPtCut)continue;
1862         list->Add(part);
1863         iCount++;
1864       }
1865       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1866         if(part->Charge()==0)continue;
1867         if(part->Pt()<fTrackPtCut)continue;
1868         if(kTrackAODMCCharged){
1869           list->Add(part);
1870         }
1871         else {
1872           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1873           list->Add(part);
1874         }
1875         iCount++;
1876       }
1877     }
1878   }// AODMCparticle
1879   list->Sort();
1880   return iCount;
1881 }
1882
1883 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1884
1885   TFile *f = TFile::Open(fPathTrPtResolution.Data());
1886
1887   if(!f)return;
1888
1889   TProfile *fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1890   TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1891   TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1892
1893   SetSmearResolution(kTRUE);
1894   SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
1895
1896
1897 }
1898
1899 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1900
1901   TFile *f = TFile::Open(fPathTrEfficiency.Data());
1902   if(!f)return;
1903
1904   TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1905   TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1906   TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1907
1908   SetDiceEfficiency(kTRUE);
1909
1910   if(fChangeEfficiencyFraction>0.) {
1911
1912     TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1913
1914     for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1915       Double_t content = hEffPosGlobSt->GetBinContent(bin);
1916       hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1917     }
1918
1919     SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1920
1921   }
1922   else
1923     SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1924
1925 }
1926
1927 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1928
1929   //
1930   // set mom res profiles
1931   //
1932
1933   if(fMomResH1) delete fMomResH1;
1934   if(fMomResH2) delete fMomResH2;
1935   if(fMomResH3) delete fMomResH3;
1936
1937   fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
1938   fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
1939   fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
1940
1941 }
1942
1943 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1944   //
1945   // set tracking efficiency histos
1946   //
1947
1948   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1949   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1950   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1951 }
1952
1953 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1954
1955   //
1956   // Get smearing on generated momentum
1957   //
1958
1959   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1960
1961   TProfile *fMomRes = 0x0;
1962   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1963   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1964   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1965
1966   if(!fMomRes) {
1967     return 0.;
1968   }
1969
1970
1971   Double_t smear = 0.;
1972
1973   if(pt>20.) {
1974     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1975     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1976     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1977   }
1978   else {
1979
1980     Int_t bin = fMomRes->FindBin(pt);
1981
1982     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1983
1984   }
1985
1986   if(fMomRes) delete fMomRes;
1987
1988   return smear;
1989 }
1990
1991 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1992   //
1993   // Fit linear function on momentum resolution at high pT
1994   //
1995
1996   if(!fMomResH1Fit && fMomResH1) {
1997     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1998     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1999     fMomResH1Fit ->SetRange(5.,100.);
2000   }
2001
2002   if(!fMomResH2Fit && fMomResH2) {
2003     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2004     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2005     fMomResH2Fit ->SetRange(5.,100.);
2006   }
2007
2008   if(!fMomResH3Fit && fMomResH3) {
2009     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2010     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2011     fMomResH3Fit ->SetRange(5.,100.);
2012   }
2013
2014 }
2015
2016 /*
2017 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2018   for(int i = 0; i < particles.GetEntries(); i++){
2019     AliVParticle *vp = (AliVParticle*)particles.At(i);
2020     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2021     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2022     jInp.set_user_index(i);
2023     inputParticles.push_back(jInp);
2024   }
2025
2026   return 0;
2027
2028 }
2029 */