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