JETAN/DEV/AliAnalysisTaskJetCluster.cxx

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