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