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