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