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   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 632   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 633   fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
 634
 635   inclusiveJets = clustSeq.inclusive_jets();
 636   sortedJets    = sorted_by_pt(inclusiveJets);
 637
 638   fh1NJetsRec->Fill(sortedJets.size());
 639
 640  // loop over all jets an fill information, first on is the leading jet
 641
 642   Int_t nRecOver = inclusiveJets.size();
 643   Int_t nRec     = inclusiveJets.size();
 644   if(inclusiveJets.size()>0){
 645     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 646     Float_t pt = leadingJet.Pt();
 647     Int_t nAodOutJets = 0;
 648     Int_t nAodOutTracks = 0;
 649     AliAODJet *aodOutJet = 0;
 650
 651     Int_t iCount = 0;
 652     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 653       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 654       while(pt<ptCut&&iCount<nRec){
 655         nRecOver--;
 656         iCount++;
 657         if(iCount<nRec){
 658           pt = sortedJets[iCount].perp();
 659         }
 660       }
 661       if(nRecOver<=0)break;
 662       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 663     }
 664     Float_t phi = leadingJet.Phi();
 665     if(phi<0)phi+=TMath::Pi()*2.;
 666     Float_t eta = leadingJet.Eta();
 667     pt = leadingJet.Pt();
 668
 669     // correlation of leading jet with tracks
 670     TIterator *recIter = recParticles.MakeIterator();
 671     recIter->Reset();
 672     AliVParticle *tmpRecTrack = 0;
 673     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 674       Float_t tmpPt = tmpRecTrack->Pt();
 675       // correlation
 676       Float_t tmpPhi =  tmpRecTrack->Phi();
 677       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 678       Float_t dPhi = phi - tmpPhi;
 679       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 680       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 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         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
 850         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
 851       }
 852
 853
 854
 855     for(int j = 0; j < nRecRan;j++){
 856       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
 857       Float_t tmpPt = tmpRec.Pt();
 858       fh1PtJetsRecInRan->Fill(tmpPt);
 859       // Fill Spectra with constituents
 860       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
 861       fh1NConstRecRan->Fill(constituents.size());
 862       fh2NConstPtRan->Fill(tmpPt,constituents.size());
 863       fh2PtNchRan->Fill(nCh,tmpPt);
 864       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
 865       // correlation
 866       Float_t tmpPhi =  tmpRec.Phi();
 867       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 868
 869       if(j==0){
 870         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
 871         fh1NConstLeadingRecRan->Fill(constituents.size());
 872         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
 873         continue;
 874       }
 875     }
 876  }
 877
 878
 879  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 880  PostData(1, fHistList);
 881 }
 882
 883 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
 884 {
 885 // Terminate analysis
 886 //
 887     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
 888 }
 889
 890
 891 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
 892
 893   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 894
 895   Int_t iCount = 0;
 896   if(type==kTrackAOD){
 897     AliAODEvent *aod = 0;
 898     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 899     else aod = AODEvent();
 900     if(!aod){
 901       return iCount;
 902     }
 903     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 904       AliAODTrack *tr = aod->GetTrack(it);
 905       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 906       if(TMath::Abs(tr->Eta())>0.9)continue;
 907       //      if(tr->Pt()<0.3)continue;
 908       list->Add(tr);
 909       iCount++;
 910     }
 911   }
 912   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 913     AliMCEvent* mcEvent = MCEvent();
 914     if(!mcEvent)return iCount;
 915     // we want to have alivpartilces so use get track
 916     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 917       if(!mcEvent->IsPhysicalPrimary(it))continue;
 918       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 919       if(type == kTrackKineAll){
 920         list->Add(part);
 921         iCount++;
 922       }
 923       else if(type == kTrackKineCharged){
 924         if(part->Particle()->GetPDG()->Charge()==0)continue;
 925         list->Add(part);
 926         iCount++;
 927       }
 928     }
 929   }
 930   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 931     AliAODEvent *aod = 0;
 932     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 933     else aod = AODEvent();
 934     if(!aod)return iCount;
 935     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 936     if(!tca)return iCount;
 937     for(int it = 0;it < tca->GetEntriesFast();++it){
 938       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
 939       if(!part->IsPhysicalPrimary())continue;
 940       if(type == kTrackAODMCAll){
 941         list->Add(part);
 942         iCount++;
 943       }
 944       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 945         if(part->Charge()==0)continue;
 946         if(kTrackAODMCCharged){
 947           list->Add(part);
 948         }
 949         else {
 950           if(TMath::Abs(part->Eta())>0.9)continue;
 951           list->Add(part);
 952         }
 953         iCount++;
 954       }
 955     }
 956   }// AODMCparticle
 957   list->Sort();
 958   return iCount;
 959
 960 }
 961
 962 /*
 963 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
 964   for(int i = 0; i < particles.GetEntries(); i++){
 965     AliVParticle *vp = (AliVParticle*)particles.At(i);
 966     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 967     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 968     jInp.set_user_index(i);
 969     inputParticles.push_back(jInp);
 970   }
 971
 972   return 0;
 973
 974 }
 975 */