JETAN/AliAnalysisTaskJetCluster.cxx

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