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