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