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   fTrackEtaWindow(0.9),
 102   fRecEtaWindow(0.5),
 103   fTrackPtCut(0.),
 104   fJetOutputMinPt(1),
 105   fJetTriggerPtCut(0),
 106   fCentCutUp(0),
 107   fCentCutLo(0),
 108   fNonStdBranch(""),
 109   fBackgroundBranch(""),
 110   fNonStdFile(""),
 111   fRparam(1.0),
 112   fAlgorithm(fastjet::kt_algorithm),
 113   fStrategy(fastjet::Best),
 114   fRecombScheme(fastjet::BIpt_scheme),
 115   fAreaType(fastjet::active_area),
 116   fGhostArea(0.01),
 117   fActiveAreaRepeats(1),
 118   fGhostEtamax(1.5),
 119   fTCAJetsOut(0x0),
 120   fTCAJetsOutRan(0x0),
 121   fTCARandomConesOut(0x0),
 122   fTCARandomConesOutRan(0x0),
 123   fAODJetBackgroundOut(0x0),
 124   fRandom(0),
 125   fh1Xsec(0x0),
 126   fh1Trials(0x0),
 127   fh1PtHard(0x0),
 128   fh1PtHardNoW(0x0),
 129   fh1PtHardTrials(0x0),
 130   fh1NJetsRec(0x0),
 131   fh1NConstRec(0x0),
 132   fh1NConstLeadingRec(0x0),
 133   fh1PtJetsRecIn(0x0),
 134   fh1PtJetsLeadingRecIn(0x0),
 135   fh1PtJetConstRec(0x0),
 136   fh1PtJetConstLeadingRec(0x0),
 137   fh1PtTracksRecIn(0x0),
 138   fh1PtTracksLeadingRecIn(0x0),
 139   fh1NJetsRecRan(0x0),
 140   fh1NConstRecRan(0x0),
 141   fh1PtJetsLeadingRecInRan(0x0),
 142   fh1NConstLeadingRecRan(0x0),
 143   fh1PtJetsRecInRan(0x0),
 144   fh1PtTracksGenIn(0x0),
 145   fh1Nch(0x0),
 146   fh1CentralityPhySel(0x0),
 147   fh1Centrality(0x0),
 148   fh1CentralitySelect(0x0),
 149   fh1ZPhySel(0x0),
 150   fh1Z(0x0),
 151   fh1ZSelect(0x0),
 152   fh2NRecJetsPt(0x0),
 153   fh2NRecTracksPt(0x0),
 154   fh2NConstPt(0x0),
 155   fh2NConstLeadingPt(0x0),
 156   fh2JetPhiEta(0x0),
 157   fh2LeadingJetPhiEta(0x0),
 158   fh2JetEtaPt(0x0),
 159   fh2LeadingJetEtaPt(0x0),
 160   fh2TrackEtaPt(0x0),
 161   fh2LeadingTrackEtaPt(0x0),
 162   fh2JetsLeadingPhiEta(0x0),
 163   fh2JetsLeadingPhiPt(0x0),
 164   fh2TracksLeadingPhiEta(0x0),
 165   fh2TracksLeadingPhiPt(0x0),
 166   fh2TracksLeadingJetPhiPt(0x0),
 167   fh2JetsLeadingPhiPtW(0x0),
 168   fh2TracksLeadingPhiPtW(0x0),
 169   fh2TracksLeadingJetPhiPtW(0x0),
 170   fh2NRecJetsPtRan(0x0),
 171   fh2NConstPtRan(0x0),
 172   fh2NConstLeadingPtRan(0x0),
 173   fh2PtNch(0x0),
 174   fh2PtNchRan(0x0),
 175   fh2PtNchN(0x0),
 176   fh2PtNchNRan(0x0),
 177   fh2TracksLeadingJetPhiPtRan(0x0),
 178   fh2TracksLeadingJetPhiPtWRan(0x0),
 179   fHistList(0x0)
 180 {
 181   for(int i = 0;i<3;i++){
 182     fh1BiARandomCones[i] = 0;
 183     fh1BiARandomConesRan[i] = 0;
 184   }
 185   for(int i = 0;i<kMaxCent;i++){
 186     fh2JetsLeadingPhiPtC[i] = 0;
 187     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 188     fh2TracksLeadingJetPhiPtC[i] = 0;
 189     fh2TracksLeadingJetPhiPtWC[i] = 0;
 190   }
 191 }
 192
 193 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
 194   AliAnalysisTaskSE(name),
 195   fAOD(0x0),
 196   fAODExtension(0x0),
 197   fRef(new TRefArray),
 198   fUseAODTrackInput(kFALSE),
 199   fUseAODMCInput(kFALSE),
 200   fUseGlobalSelection(kFALSE),
 201   fUseBackgroundCalc(kFALSE),
 202   fFilterMask(0),
 203   fTrackTypeRec(kTrackUndef),
 204   fTrackTypeGen(kTrackUndef),
 205   fNSkipLeadingRan(0),
 206   fNRandomCones(0),
 207   fAvgTrials(1),
 208   fExternalWeight(1),
 209   fTrackEtaWindow(0.9),
 210   fRecEtaWindow(0.5),
 211   fTrackPtCut(0.),
 212   fJetOutputMinPt(1),
 213   fJetTriggerPtCut(0),
 214   fCentCutUp(0),
 215   fCentCutLo(0),
 216   fNonStdBranch(""),
 217   fBackgroundBranch(""),
 218   fNonStdFile(""),
 219   fRparam(1.0),
 220   fAlgorithm(fastjet::kt_algorithm),
 221   fStrategy(fastjet::Best),
 222   fRecombScheme(fastjet::BIpt_scheme),
 223   fAreaType(fastjet::active_area),
 224   fGhostArea(0.01),
 225   fActiveAreaRepeats(1),
 226   fGhostEtamax(1.5),
 227   fTCAJetsOut(0x0),
 228   fTCAJetsOutRan(0x0),
 229   fTCARandomConesOut(0x0),
 230   fTCARandomConesOutRan(0x0),
 231   fAODJetBackgroundOut(0x0),
 232   fRandom(0),
 233   fh1Xsec(0x0),
 234   fh1Trials(0x0),
 235   fh1PtHard(0x0),
 236   fh1PtHardNoW(0x0),
 237   fh1PtHardTrials(0x0),
 238   fh1NJetsRec(0x0),
 239   fh1NConstRec(0x0),
 240   fh1NConstLeadingRec(0x0),
 241   fh1PtJetsRecIn(0x0),
 242   fh1PtJetsLeadingRecIn(0x0),
 243   fh1PtJetConstRec(0x0),
 244   fh1PtJetConstLeadingRec(0x0),
 245   fh1PtTracksRecIn(0x0),
 246   fh1PtTracksLeadingRecIn(0x0),
 247   fh1NJetsRecRan(0x0),
 248   fh1NConstRecRan(0x0),
 249   fh1PtJetsLeadingRecInRan(0x0),
 250   fh1NConstLeadingRecRan(0x0),
 251   fh1PtJetsRecInRan(0x0),
 252   fh1PtTracksGenIn(0x0),
 253   fh1Nch(0x0),
 254   fh1CentralityPhySel(0x0),
 255   fh1Centrality(0x0),
 256   fh1CentralitySelect(0x0),
 257   fh1ZPhySel(0x0),
 258   fh1Z(0x0),
 259   fh1ZSelect(0x0),
 260   fh2NRecJetsPt(0x0),
 261   fh2NRecTracksPt(0x0),
 262   fh2NConstPt(0x0),
 263   fh2NConstLeadingPt(0x0),
 264   fh2JetPhiEta(0x0),
 265   fh2LeadingJetPhiEta(0x0),
 266   fh2JetEtaPt(0x0),
 267   fh2LeadingJetEtaPt(0x0),
 268   fh2TrackEtaPt(0x0),
 269   fh2LeadingTrackEtaPt(0x0),
 270   fh2JetsLeadingPhiEta(0x0),
 271   fh2JetsLeadingPhiPt(0x0),
 272   fh2TracksLeadingPhiEta(0x0),
 273   fh2TracksLeadingPhiPt(0x0),
 274   fh2TracksLeadingJetPhiPt(0x0),
 275   fh2JetsLeadingPhiPtW(0x0),
 276   fh2TracksLeadingPhiPtW(0x0),
 277   fh2TracksLeadingJetPhiPtW(0x0),
 278   fh2NRecJetsPtRan(0x0),
 279   fh2NConstPtRan(0x0),
 280   fh2NConstLeadingPtRan(0x0),
 281   fh2PtNch(0x0),
 282   fh2PtNchRan(0x0),
 283   fh2PtNchN(0x0),
 284   fh2PtNchNRan(0x0),
 285   fh2TracksLeadingJetPhiPtRan(0x0),
 286   fh2TracksLeadingJetPhiPtWRan(0x0),
 287   fHistList(0x0)
 288 {
 289   for(int i = 0;i<3;i++){
 290     fh1BiARandomCones[i] = 0;
 291     fh1BiARandomConesRan[i] = 0;
 292   }
 293   for(int i = 0;i<kMaxCent;i++){
 294     fh2JetsLeadingPhiPtC[i] = 0;
 295     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 296     fh2TracksLeadingJetPhiPtC[i] = 0;
 297     fh2TracksLeadingJetPhiPtWC[i] = 0;
 298   }
 299  DefineOutput(1, TList::Class());
 300 }
 301
 302
 303
 304 Bool_t AliAnalysisTaskJetCluster::Notify()
 305 {
 306   //
 307   // Implemented Notify() to read the cross sections
 308   // and number of trials from pyxsec.root
 309   //
 310   return kTRUE;
 311 }
 312
 313 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 314 {
 315
 316   //
 317   // Create the output container
 318   //
 319
 320   fRandom = new TRandom3(0);
 321
 322
 323   // Connect the AOD
 324
 325
 326   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 327
 328
 329
 330   if(fNonStdBranch.Length()!=0)
 331     {
 332       // only create the output branch if we have a name
 333       // Create a new branch for jets...
 334       //  -> cleared in the UserExec....
 335       // here we can also have the case that the brnaches are written to a separate file
 336
 337       fTCAJetsOut = new TClonesArray("AliAODJet", 0);
 338       fTCAJetsOut->SetName(fNonStdBranch.Data());
 339       AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
 340
 341
 342       fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
 343       fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
 344       AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
 345
 346       if(fUseBackgroundCalc){
 347         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
 348           fAODJetBackgroundOut = new AliAODJetEventBackground();
 349           fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
 350           AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
 351         }
 352       }
 353       // create the branch for the random cones with the same R
 354       TString cName = Form("%sRandomCone",fNonStdBranch.Data());
 355
 356       if(fNRandomCones>0){
 357         if(!AODEvent()->FindListObject(cName.Data())){
 358           fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
 359           fTCARandomConesOut->SetName(cName.Data());
 360           AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
 361         }
 362         // create the branch with the random for the random cones on the random event
 363         cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 364         if(!AODEvent()->FindListObject(cName.Data())){
 365           fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
 366           fTCARandomConesOutRan->SetName(cName.Data());
 367           AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
 368         }
 369       }
 370
 371       if(fNonStdFile.Length()!=0){
 372         //
 373         // case that we have an AOD extension we need to fetch the jets from the extended output
 374         // we identify the extension aod event by looking for the branchname
 375         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 376         // case that we have an AOD extension we need can fetch the background maybe from the extended output
 377         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 378       }
 379     }
 380
 381
 382   if(!fHistList)fHistList = new TList();
 383   fHistList->SetOwner();
 384
 385   Bool_t oldStatus = TH1::AddDirectoryStatus();
 386   TH1::AddDirectory(kFALSE);
 387
 388   //
 389   //  Histogram
 390
 391   const Int_t nBinPt = 200;
 392   Double_t binLimitsPt[nBinPt+1];
 393   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 394     if(iPt == 0){
 395       binLimitsPt[iPt] = 0.0;
 396     }
 397     else {// 1.0
 398       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 399     }
 400   }
 401
 402   const Int_t nBinPhi = 90;
 403   Double_t binLimitsPhi[nBinPhi+1];
 404   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 405     if(iPhi==0){
 406       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 407     }
 408     else{
 409       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 410     }
 411   }
 412
 413
 414
 415   const Int_t nBinEta = 40;
 416   Double_t binLimitsEta[nBinEta+1];
 417   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 418     if(iEta==0){
 419       binLimitsEta[iEta] = -2.0;
 420     }
 421     else{
 422       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 423     }
 424   }
 425
 426   const int nChMax = 4000;
 427
 428   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 429   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 430
 431   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 432   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 433
 434
 435   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 436   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 437
 438   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 439   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 440   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 441   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 442
 443
 444   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 445   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 446   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 447
 448   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 449   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 450   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 451   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 452   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 453   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 454   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 455   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 456   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 457   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 458
 459   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 460   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
 461   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 462
 463   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 464   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
 465   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 466
 467   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 468   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 469   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 470   //
 471
 472
 473   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 474   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 475   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 476   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 477
 478   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 479   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);
 480   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);
 481   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);
 482
 483
 484
 485   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 486                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 487   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 488                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 489
 490   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 491                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 492   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 493                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 494
 495   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 496                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 497   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 498                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 499
 500
 501
 502   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 503                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 504   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 505                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 506   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 507                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 508   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 509                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 510   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 511                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 512   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 513                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 514
 515   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 516                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 517   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 518                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 519
 520   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 521                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 522   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 523                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 524
 525
 526   if(fNRandomCones>0&&fUseBackgroundCalc){
 527     for(int i = 0;i<3;i++){
 528       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 529       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 530     }
 531   }
 532
 533   for(int i = 0;i < kMaxCent;i++){
 534     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
 535     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
 536     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
 537     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
 538   }
 539
 540   const Int_t saveLevel = 3; // large save level more histos
 541   if(saveLevel>0){
 542     fHistList->Add(fh1Xsec);
 543     fHistList->Add(fh1Trials);
 544
 545     fHistList->Add(fh1NJetsRec);
 546     fHistList->Add(fh1NConstRec);
 547     fHistList->Add(fh1NConstLeadingRec);
 548     fHistList->Add(fh1PtJetsRecIn);
 549     fHistList->Add(fh1PtJetsLeadingRecIn);
 550     fHistList->Add(fh1PtTracksRecIn);
 551     fHistList->Add(fh1PtTracksLeadingRecIn);
 552     fHistList->Add(fh1PtJetConstRec);
 553     fHistList->Add(fh1PtJetConstLeadingRec);
 554     fHistList->Add(fh1NJetsRecRan);
 555     fHistList->Add(fh1NConstRecRan);
 556     fHistList->Add(fh1PtJetsLeadingRecInRan);
 557     fHistList->Add(fh1NConstLeadingRecRan);
 558     fHistList->Add(fh1PtJetsRecInRan);
 559     fHistList->Add(fh1Nch);
 560     fHistList->Add(fh1Centrality);
 561     fHistList->Add(fh1CentralitySelect);
 562     fHistList->Add(fh1CentralityPhySel);
 563     fHistList->Add(fh1Z);
 564     fHistList->Add(fh1ZSelect);
 565     fHistList->Add(fh1ZPhySel);
 566     if(fNRandomCones>0&&fUseBackgroundCalc){
 567       for(int i = 0;i<3;i++){
 568         fHistList->Add(fh1BiARandomCones[i]);
 569         fHistList->Add(fh1BiARandomConesRan[i]);
 570       }
 571     }
 572   for(int i = 0;i < kMaxCent;i++){
 573     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
 574     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
 575     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
 576     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
 577   }
 578
 579     fHistList->Add(fh2NRecJetsPt);
 580     fHistList->Add(fh2NRecTracksPt);
 581     fHistList->Add(fh2NConstPt);
 582     fHistList->Add(fh2NConstLeadingPt);
 583     fHistList->Add(fh2PtNch);
 584     fHistList->Add(fh2PtNchRan);
 585     fHistList->Add(fh2PtNchN);
 586     fHistList->Add(fh2PtNchNRan);
 587     fHistList->Add(fh2JetPhiEta);
 588     fHistList->Add(fh2LeadingJetPhiEta);
 589     fHistList->Add(fh2JetEtaPt);
 590     fHistList->Add(fh2LeadingJetEtaPt);
 591     fHistList->Add(fh2TrackEtaPt);
 592     fHistList->Add(fh2LeadingTrackEtaPt);
 593     fHistList->Add(fh2JetsLeadingPhiEta );
 594     fHistList->Add(fh2JetsLeadingPhiPt);
 595     fHistList->Add(fh2TracksLeadingPhiEta);
 596     fHistList->Add(fh2TracksLeadingPhiPt);
 597     fHistList->Add(fh2TracksLeadingJetPhiPt);
 598     fHistList->Add(fh2JetsLeadingPhiPtW);
 599     fHistList->Add(fh2TracksLeadingPhiPtW);
 600     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 601     fHistList->Add(fh2NRecJetsPtRan);
 602     fHistList->Add(fh2NConstPtRan);
 603     fHistList->Add(fh2NConstLeadingPtRan);
 604     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 605     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 606   }
 607
 608   // =========== Switch on Sumw2 for all histos ===========
 609   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 610     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 611     if (h1){
 612       h1->Sumw2();
 613       continue;
 614     }
 615     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 616     if(hn)hn->Sumw2();
 617   }
 618   TH1::AddDirectory(oldStatus);
 619 }
 620
 621 void AliAnalysisTaskJetCluster::Init()
 622 {
 623   //
 624   // Initialization
 625   //
 626
 627   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 628
 629 }
 630
 631 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 632 {
 633
 634   if(fUseGlobalSelection){
 635     // no selection by the service task, we continue
 636     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 637     PostData(1, fHistList);
 638     return;
 639   }
 640
 641
 642
 643   // handle and reset the output jet branch
 644
 645   if(fTCAJetsOut)fTCAJetsOut->Delete();
 646   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
 647   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
 648   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
 649   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
 650
 651   AliAODJetEventBackground* externalBackground = 0;
 652   if(!externalBackground&&fBackgroundBranch.Length()){
 653     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
 654     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 655   }
 656   //
 657   // Execute analysis for current event
 658   //
 659   AliESDEvent *fESD = 0;
 660   if(fUseAODTrackInput){
 661     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 662     if(!fAOD){
 663       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 664       return;
 665     }
 666     // fethc the header
 667   }
 668   else{
 669     //  assume that the AOD is in the general output...
 670     fAOD  = AODEvent();
 671     if(!fAOD){
 672       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 673       return;
 674     }
 675     if(fDebug>0){
 676       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 677     }
 678   }
 679
 680   Bool_t selectEvent =  false;
 681   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 682
 683   Float_t cent = 0;
 684   Float_t zVtx  = 0;
 685   Int_t cenClass = -1;
 686   if(fAOD){
 687     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 688     TString vtxTitle(vtxAOD->GetTitle());
 689     zVtx = vtxAOD->GetZ();
 690
 691     cent = fAOD->GetHeader()->GetCentrality();
 692     if(cent<10)cenClass = 0;
 693     else if(cent<30)cenClass = 1;
 694     else if(cent<50)cenClass = 2;
 695     else if(cent<80)cenClass = 3;
 696     if(physicsSelection){
 697       fh1CentralityPhySel->Fill(cent);
 698       fh1ZPhySel->Fill(zVtx);
 699     }
 700
 701
 702     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 703         Float_t yvtx = vtxAOD->GetY();
 704         Float_t xvtx = vtxAOD->GetX();
 705         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 706         if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
 707           if(physicsSelection){
 708             selectEvent = true;
 709           }
 710         }
 711     }
 712     if(fCentCutUp>0){
 713       if(cent<fCentCutLo||cent>fCentCutUp){
 714         selectEvent = false;
 715       }
 716     }
 717
 718   }
 719   if(!selectEvent){
 720     PostData(1, fHistList);
 721     return;
 722   }
 723   fh1Centrality->Fill(cent);
 724   fh1Z->Fill(zVtx);
 725   fh1Trials->Fill("#sum{ntrials}",1);
 726
 727
 728   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 729
 730   // ==== General variables needed
 731
 732
 733
 734   // we simply fetch the tracks/mc particles as a list of AliVParticles
 735
 736   TList recParticles;
 737   TList genParticles;
 738
 739   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 740   Float_t nCh = recParticles.GetEntries();
 741   fh1Nch->Fill(nCh);
 742   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 743   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 744   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 745
 746   // find the jets....
 747
 748   vector<fastjet::PseudoJet> inputParticlesRec;
 749   vector<fastjet::PseudoJet> inputParticlesRecRan;
 750
 751   // Generate the random cones
 752
 753   AliAODJet vTmpRan(1,0,0,1);
 754   for(int i = 0; i < recParticles.GetEntries(); i++){
 755     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 756     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 757     // we take total momentum here
 758     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 759     jInp.set_user_index(i);
 760     inputParticlesRec.push_back(jInp);
 761
 762     // the randomized input changes eta and phi, but keeps the p_T
 763     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 764       Double_t pT = vp->Pt();
 765       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
 766       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
 767
 768       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
 769       Double_t pZ = pT/TMath::Tan(theta);
 770
 771       Double_t pX = pT * TMath::Cos(phi);
 772       Double_t pY = pT * TMath::Sin(phi);
 773       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 774       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 775
 776       jInpRan.set_user_index(i);
 777       inputParticlesRecRan.push_back(jInpRan);
 778       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
 779     }
 780
 781     // fill the tref array, only needed when we write out jets
 782     if(fTCAJetsOut){
 783       if(i == 0){
 784         fRef->Delete(); // make sure to delete before placement new...
 785         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
 786       }
 787       fRef->Add(vp);
 788     }
 789   }// recparticles
 790
 791   if(inputParticlesRec.size()==0){
 792     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 793     PostData(1, fHistList);
 794     return;
 795   }
 796
 797   // run fast jet
 798   // employ setters for these...
 799
 800
 801   // now create the object that holds info about ghosts
 802   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
 803     // reduce CPU time...
 804     fGhostArea = 0.5;
 805     fActiveAreaRepeats = 0;
 806   }
 807
 808  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 809   fastjet::AreaType areaType =   fastjet::active_area;
 810   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 811   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 812   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 813   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
 814
 815   //range where to compute background
 816   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 817   phiMin = 0;
 818   phiMax = 2*TMath::Pi();
 819   rapMax = fGhostEtamax - fRparam;
 820   rapMin = - fGhostEtamax + fRparam;
 821   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 822
 823
 824
 825   inclusiveJets = clustSeq.inclusive_jets();
 826   sortedJets    = sorted_by_pt(inclusiveJets);
 827
 828   fh1NJetsRec->Fill(sortedJets.size());
 829
 830  // loop over all jets an fill information, first on is the leading jet
 831
 832   Int_t nRecOver = inclusiveJets.size();
 833   Int_t nRec     = inclusiveJets.size();
 834   if(inclusiveJets.size()>0){
 835     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 836     Double_t area = clustSeq.area(sortedJets[0]);
 837     leadingJet.SetEffArea(area,0);
 838     Float_t pt = leadingJet.Pt();
 839     Int_t nAodOutJets = 0;
 840     Int_t nAodOutTracks = 0;
 841     AliAODJet *aodOutJet = 0;
 842
 843     Int_t iCount = 0;
 844     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 845       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 846       while(pt<ptCut&&iCount<nRec){
 847         nRecOver--;
 848         iCount++;
 849         if(iCount<nRec){
 850           pt = sortedJets[iCount].perp();
 851         }
 852       }
 853       if(nRecOver<=0)break;
 854       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 855     }
 856     Float_t phi = leadingJet.Phi();
 857     if(phi<0)phi+=TMath::Pi()*2.;
 858     Float_t eta = leadingJet.Eta();
 859     Float_t pTback = 0;
 860     if(externalBackground){
 861       // carefull has to be filled in a task before
 862       // todo, ReArrange to the botom
 863       pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
 864     }
 865     pt = leadingJet.Pt() - pTback;
 866     // correlation of leading jet with tracks
 867     TIterator *recIter = recParticles.MakeIterator();
 868     recIter->Reset();
 869     AliVParticle *tmpRecTrack = 0;
 870     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 871       Float_t tmpPt = tmpRecTrack->Pt();
 872       // correlation
 873       Float_t tmpPhi =  tmpRecTrack->Phi();
 874       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 875       Float_t dPhi = phi - tmpPhi;
 876       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 877       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 878       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 879       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 880       if(cenClass>=0){
 881         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
 882         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
 883       }
 884
 885     }
 886
 887
 888     TLorentzVector vecareab;
 889     for(int j = 0; j < nRec;j++){
 890       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 891       aodOutJet = 0;
 892       nAodOutTracks = 0;
 893       Float_t tmpPt = tmpRec.Pt();
 894
 895       if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
 896         aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
 897         Double_t area1 = clustSeq.area(sortedJets[j]);
 898         aodOutJet->SetEffArea(area1,0);
 899         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
 900         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
 901         aodOutJet->SetVectorAreaCharged(&vecareab);
 902
 903
 904       }
 905
 906
 907       Float_t tmpPtBack = 0;
 908       if(externalBackground){
 909         // carefull has to be filled in a task before
 910        // todo, ReArrange to the botom
 911         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
 912       }
 913       tmpPt = tmpPt - tmpPtBack;
 914       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
 915
 916       fh1PtJetsRecIn->Fill(tmpPt);
 917       // Fill Spectra with constituents
 918       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 919
 920       fh1NConstRec->Fill(constituents.size());
 921       fh2PtNch->Fill(nCh,tmpPt);
 922       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 923       fh2NConstPt->Fill(tmpPt,constituents.size());
 924       // loop over constiutents and fill spectrum
 925
 926      for(unsigned int ic = 0; ic < constituents.size();ic++){
 927        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 928        fh1PtJetConstRec->Fill(part->Pt());
 929        if(aodOutJet){
 930          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
 931        }
 932        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
 933      }
 934
 935      // correlation
 936      Float_t tmpPhi =  tmpRec.Phi();
 937      Float_t tmpEta =  tmpRec.Eta();
 938      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 939      if(j==0){
 940        fh1PtJetsLeadingRecIn->Fill(tmpPt);
 941        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
 942        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
 943        fh1NConstLeadingRec->Fill(constituents.size());
 944        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
 945        continue;
 946      }
 947      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
 948      fh2JetEtaPt->Fill(tmpEta,tmpPt);
 949      Float_t dPhi = phi - tmpPhi;
 950      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 951      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 952      Float_t dEta = eta - tmpRec.Eta();
 953      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 954      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 955      if(cenClass>=0){
 956        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
 957        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
 958      }
 959      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 960     }// loop over reconstructed jets
 961    delete recIter;
 962
 963
 964
 965    // Add the random cones...
 966    if(fNRandomCones>0&&fTCARandomConesOut){
 967      // create a random jet within the acceptance
 968      Double_t etaMax = 0.8 - fRparam;
 969      Int_t nCone = 0;
 970      Int_t nConeRan = 0;
 971      Double_t pTC = 1; // small number
 972      for(int ir = 0;ir < fNRandomCones;ir++){
 973        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
 974        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
 975        // massless jet
 976        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 977        Double_t pZC = pTC/TMath::Tan(thetaC);
 978        Double_t pXC = pTC * TMath::Cos(phiC);
 979        Double_t pYC = pTC * TMath::Sin(phiC);
 980        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 981        AliAODJet tmpRecC (pXC,pYC,pZC, pC);
 982        bool skip = false;
 983        for(int jj = 0; jj < TMath::Min(nRec,2);jj++){// test for overlap with leading jets
 984          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
 985          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
 986            skip = true;
 987            break;
 988          }
 989        }
 990        // test for overlap with previous cones to avoid double counting
 991        for(int iic = 0;iic<ir;iic++){
 992          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
 993          if(iicone){
 994            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
 995              skip = true;
 996              break;
 997            }
 998          }
 999        }
1000        if(skip)continue;
1001        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1002        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nConeRan++]) AliAODJet(tmpRecC);
1003        if(fTCARandomConesOutRan)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 = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
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())>fTrackEtaWindow)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())>fTrackEtaWindow)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 */