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