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     AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
 672
 673     recIter->Reset();
 674     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 675       Float_t tmpPt = tmpRecTrack->Pt();
 676       // correlation
 677       Float_t tmpPhi =  tmpRecTrack->Phi();
 678       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 679       Float_t dPhi = phi - tmpPhi;
 680       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 681       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 682       //      Float_t dEta = eta - tmpRecTrack->Eta();
 683       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 684       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 685     }
 686
 687
 688
 689    for(int j = 0; j < nRec;j++){
 690      AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 691      aodOutJet = 0;
 692      nAodOutTracks = 0;
 693      Float_t tmpPt = tmpRec.Pt();
 694      fh1PtJetsRecIn->Fill(tmpPt);
 695      // Fill Spectra with constituents
 696      vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 697      fh1NConstRec->Fill(constituents.size());
 698      fh2PtNch->Fill(nCh,tmpPt);
 699      fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 700      fh2NConstPt->Fill(tmpPt,constituents.size());
 701      // loop over constiutents and fill spectrum
 702
 703      // Add the jet information and the track references to the output AOD
 704
 705      if(tmpPt>fJetOutputMinPt&&jarray){
 706        aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
 707      }
 708
 709
 710      for(unsigned int ic = 0; ic < constituents.size();ic++){
 711        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
 712        fh1PtJetConstRec->Fill(part->Pt());
 713        if(aodOutJet){
 714          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
 715        }
 716        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
 717      }
 718
 719
 720
 721
 722
 723
 724
 725
 726      // correlation
 727      Float_t tmpPhi =  tmpRec.Phi();
 728      Float_t tmpEta =  tmpRec.Eta();
 729      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 730
 731      if(j==0){
 732        fh1PtJetsLeadingRecIn->Fill(tmpPt);
 733        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
 734        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
 735        fh1NConstLeadingRec->Fill(constituents.size());
 736        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
 737        continue;
 738      }
 739      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
 740      fh2JetEtaPt->Fill(tmpEta,tmpPt);
 741      Float_t dPhi = phi - tmpPhi;
 742      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 743      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 744      Float_t dEta = eta - tmpRec.Eta();
 745      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 746      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 747      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 748    }
 749    delete recIter;
 750   }
 751
 752  // fill track information
 753  Int_t nTrackOver = recParticles.GetSize();
 754   // do the same for tracks and jets
 755  if(nTrackOver>0){
 756    TIterator *recIter = recParticles.MakeIterator();
 757    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 758    Float_t pt = tmpRec->Pt();
 759    //    Printf("Leading track p_t %3.3E",pt);
 760    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 761      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 762      while(pt<ptCut&&tmpRec){
 763        nTrackOver--;
 764        tmpRec = (AliVParticle*)(recIter->Next());
 765        if(tmpRec){
 766          pt = tmpRec->Pt();
 767        }
 768      }
 769      if(nTrackOver<=0)break;
 770      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 771    }
 772
 773    recIter->Reset();
 774    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 775    Float_t phi = leading->Phi();
 776    if(phi<0)phi+=TMath::Pi()*2.;
 777    Float_t eta = leading->Eta();
 778    pt  = leading->Pt();
 779    while((tmpRec = (AliVParticle*)(recIter->Next()))){
 780      Float_t tmpPt = tmpRec->Pt();
 781      Float_t tmpEta = tmpRec->Eta();
 782      fh1PtTracksRecIn->Fill(tmpPt);
 783      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
 784      if(tmpRec==leading){
 785        fh1PtTracksLeadingRecIn->Fill(tmpPt);
 786        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
 787        continue;
 788      }
 789       // correlation
 790      Float_t tmpPhi =  tmpRec->Phi();
 791
 792      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 793      Float_t dPhi = phi - tmpPhi;
 794      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 795      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 796      Float_t dEta = eta - tmpRec->Eta();
 797      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 798      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 799      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 800    }
 801    delete recIter;
 802  }
 803
 804  // find the random jets
 805  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
 806  fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
 807
 808  inclusiveJetsRan = clustSeqRan.inclusive_jets();
 809  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
 810
 811  // fill the jet information from random track
 812
 813   fh1NJetsRecRan->Fill(sortedJetsRan.size());
 814
 815  // loop over all jets an fill information, first on is the leading jet
 816
 817  Int_t nRecOverRan = inclusiveJetsRan.size();
 818  Int_t nRecRan     = inclusiveJetsRan.size();
 819  if(inclusiveJetsRan.size()>0){
 820    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
 821    Float_t pt = leadingJet.Pt();
 822
 823    Int_t iCount = 0;
 824    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
 825      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
 826       while(pt<ptCut&&iCount<nRecRan){
 827         nRecOverRan--;
 828         iCount++;
 829         if(iCount<nRecRan){
 830           pt = sortedJetsRan[iCount].perp();
 831         }
 832       }
 833       if(nRecOverRan<=0)break;
 834       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
 835     }
 836     Float_t phi = leadingJet.Phi();
 837     if(phi<0)phi+=TMath::Pi()*2.;
 838     pt = leadingJet.Pt();
 839
 840     // correlation of leading jet with random tracks
 841
 842     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
 843       {
 844         Float_t tmpPt = inputParticlesRecRan[ip].perp();
 845         // correlation
 846         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
 847         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 848         Float_t dPhi = phi - tmpPhi;
 849         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 850         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 851         //      Float_t dEta = eta - tmpRecTrack->Eta();
 852         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
 853         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
 854       }
 855
 856
 857
 858     for(int j = 0; j < nRecRan;j++){
 859       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
 860       Float_t tmpPt = tmpRec.Pt();
 861       fh1PtJetsRecInRan->Fill(tmpPt);
 862       // Fill Spectra with constituents
 863       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
 864       fh1NConstRecRan->Fill(constituents.size());
 865       fh2NConstPtRan->Fill(tmpPt,constituents.size());
 866       fh2PtNchRan->Fill(nCh,tmpPt);
 867       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
 868       // correlation
 869       Float_t tmpPhi =  tmpRec.Phi();
 870       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 871
 872       if(j==0){
 873         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
 874         fh1NConstLeadingRecRan->Fill(constituents.size());
 875         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
 876         continue;
 877       }
 878     }
 879  }
 880
 881
 882  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 883  PostData(1, fHistList);
 884 }
 885
 886 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
 887 {
 888 // Terminate analysis
 889 //
 890     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
 891 }
 892
 893
 894 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
 895
 896   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 897
 898   Int_t iCount = 0;
 899   if(type==kTrackAOD){
 900     AliAODEvent *aod = 0;
 901     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 902     else aod = AODEvent();
 903     if(!aod){
 904       return iCount;
 905     }
 906     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 907       AliAODTrack *tr = aod->GetTrack(it);
 908       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 909       if(TMath::Abs(tr->Eta())>0.9)continue;
 910       //      if(tr->Pt()<0.3)continue;
 911       list->Add(tr);
 912       iCount++;
 913     }
 914   }
 915   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 916     AliMCEvent* mcEvent = MCEvent();
 917     if(!mcEvent)return iCount;
 918     // we want to have alivpartilces so use get track
 919     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 920       if(!mcEvent->IsPhysicalPrimary(it))continue;
 921       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 922       if(type == kTrackKineAll){
 923         list->Add(part);
 924         iCount++;
 925       }
 926       else if(type == kTrackKineCharged){
 927         if(part->Particle()->GetPDG()->Charge()==0)continue;
 928         list->Add(part);
 929         iCount++;
 930       }
 931     }
 932   }
 933   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 934     AliAODEvent *aod = 0;
 935     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 936     else aod = AODEvent();
 937     if(!aod)return iCount;
 938     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 939     if(!tca)return iCount;
 940     for(int it = 0;it < tca->GetEntriesFast();++it){
 941       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 942       if(!part->IsPhysicalPrimary())continue;
 943       if(type == kTrackAODMCAll){
 944         list->Add(part);
 945         iCount++;
 946       }
 947       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 948         if(part->Charge()==0)continue;
 949         if(kTrackAODMCCharged){
 950           list->Add(part);
 951         }
 952         else {
 953           if(TMath::Abs(part->Eta())>0.9)continue;
 954           list->Add(part);
 955         }
 956         iCount++;
 957       }
 958     }
 959   }// AODMCparticle
 960   list->Sort();
 961   return iCount;
 962
 963 }
 964
 965 /*
 966 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
 967   for(int i = 0; i < particles.GetEntries(); i++){
 968     AliVParticle *vp = (AliVParticle*)particles.At(i);
 969     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 970     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 971     jInp.set_user_index(i);
 972     inputParticles.push_back(jInp);
 973   }
 974
 975   return 0;
 976
 977 }
 978 */