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