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