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