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