JETAN/DEV/AliAnalysisTaskJetCluster.cxx

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