PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.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 <TFile.h>
  29 #include <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TList.h>
  35 #include <TLorentzVector.h>
  36 #include <TClonesArray.h>
  37 #include  "TDatabasePDG.h"
  38
  39 #include "AliAnalysisTaskJetSpectrum2.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliJetFinder.h"
  42 #include "AliJetHeader.h"
  43 #include "AliJetReader.h"
  44 #include "AliJetReaderHeader.h"
  45 #include "AliUA1JetHeaderV1.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODTrack.h"
  50 #include "AliAODJet.h"
  51 #include "AliAODMCParticle.h"
  52 #include "AliMCEventHandler.h"
  53 #include "AliMCEvent.h"
  54 #include "AliStack.h"
  55 #include "AliGenPythiaEventHeader.h"
  56 #include "AliJetKineReaderHeader.h"
  57 #include "AliGenCocktailEventHeader.h"
  58 #include "AliInputEventHandler.h"
  59
  60
  61 #include "AliAnalysisHelperJetTasks.h"
  62
  63 ClassImp(AliAnalysisTaskJetSpectrum2)
  64
  65 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
  66   fJetHeaderRec(0x0),
  67   fJetHeaderGen(0x0),
  68   fAOD(0x0),
  69   fhnCorrelation(0x0),
  70   fBranchRec("jets"),
  71   fBranchGen(""),
  72   fUseAODJetInput(kFALSE),
  73   fUseAODTrackInput(kFALSE),
  74   fUseAODMCInput(kFALSE),
  75   fUseGlobalSelection(kFALSE),
  76   fUseExternalWeightOnly(kFALSE),
  77   fLimitGenJetEta(kFALSE),
  78   fFilterMask(0),
  79   fAnalysisType(0),
  80   fTrackTypeRec(kTrackUndef),
  81   fTrackTypeGen(kTrackUndef),
  82   fAvgTrials(1),
  83   fExternalWeight(1),
  84   fRecEtaWindow(0.5),
  85   fh1Xsec(0x0),
  86   fh1Trials(0x0),
  87   fh1PtHard(0x0),
  88   fh1PtHardNoW(0x0),
  89   fh1PtHardTrials(0x0),
  90   fh1NGenJets(0x0),
  91   fh1NRecJets(0x0),
  92   fh1PtTrackRec(0x0),
  93   fh1SumPtTrackRec(0x0),
  94   fh1SumPtTrackAreaRec(0x0),
  95   fh1PtJetsRecIn(0x0),
  96   fh1PtJetsLeadingRecIn(0x0),
  97   fh1PtTracksRecIn(0x0),
  98   fh1PtTracksLeadingRecIn(0x0),
  99   fh1PtTracksGenIn(0x0),
 100   fh2NRecJetsPt(0x0),
 101   fh2NRecTracksPt(0x0),
 102   fh2JetsLeadingPhiEta(0x0),
 103   fh2JetsLeadingPhiPt(0x0),
 104   fh2TracksLeadingPhiEta(0x0),
 105   fh2TracksLeadingPhiPt(0x0),
 106   fHistList(0x0)
 107 {
 108   for(int i = 0;i < kMaxStep*2;++i){
 109     fhnJetContainer[i] = 0;
 110   }
 111   for(int i = 0;i < kMaxJets;++i){
 112     fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
 113     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 114   }
 115
 116 }
 117
 118 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 119   AliAnalysisTaskSE(name),
 120   fJetHeaderRec(0x0),
 121   fJetHeaderGen(0x0),
 122   fAOD(0x0),
 123   fhnCorrelation(0x0),
 124   fBranchRec("jets"),
 125   fBranchGen(""),
 126   fUseAODJetInput(kFALSE),
 127   fUseAODTrackInput(kFALSE),
 128   fUseAODMCInput(kFALSE),
 129   fUseGlobalSelection(kFALSE),
 130   fUseExternalWeightOnly(kFALSE),
 131   fLimitGenJetEta(kFALSE),
 132   fFilterMask(0),
 133   fAnalysisType(0),
 134   fTrackTypeRec(kTrackUndef),
 135   fTrackTypeGen(kTrackUndef),
 136   fAvgTrials(1),
 137   fExternalWeight(1),
 138   fRecEtaWindow(0.5),
 139   fh1Xsec(0x0),
 140   fh1Trials(0x0),
 141   fh1PtHard(0x0),
 142   fh1PtHardNoW(0x0),
 143   fh1PtHardTrials(0x0),
 144   fh1NGenJets(0x0),
 145   fh1NRecJets(0x0),
 146   fh1PtTrackRec(0x0),
 147   fh1SumPtTrackRec(0x0),
 148   fh1SumPtTrackAreaRec(0x0),
 149   fh1PtJetsRecIn(0x0),
 150   fh1PtJetsLeadingRecIn(0x0),
 151   fh1PtTracksRecIn(0x0),
 152   fh1PtTracksLeadingRecIn(0x0),
 153   fh1PtTracksGenIn(0x0),
 154   fh2NRecJetsPt(0x0),
 155   fh2NRecTracksPt(0x0),
 156   fh2JetsLeadingPhiEta(0x0),
 157   fh2JetsLeadingPhiPt(0x0),
 158   fh2TracksLeadingPhiEta(0x0),
 159   fh2TracksLeadingPhiPt(0x0),
 160   fHistList(0x0)
 161 {
 162
 163   for(int i = 0;i < kMaxStep*2;++i){
 164     fhnJetContainer[i] = 0;
 165   }
 166   for(int i = 0;i < kMaxJets;++i){
 167     fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
 168     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 169   }
 170   DefineOutput(1, TList::Class());
 171 }
 172
 173
 174
 175 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 176 {
 177   //
 178   // Implemented Notify() to read the cross sections
 179   // and number of trials from pyxsec.root
 180   //
 181
 182   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 183   Float_t xsection = 0;
 184   Float_t ftrials  = 1;
 185
 186   fAvgTrials = 1;
 187   if(tree){
 188     TFile *curfile = tree->GetCurrentFile();
 189     if (!curfile) {
 190       Error("Notify","No current file");
 191       return kFALSE;
 192     }
 193     if(!fh1Xsec||!fh1Trials){
 194       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 195       return kFALSE;
 196     }
 197     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 198     fh1Xsec->Fill("<#sigma>",xsection);
 199     // construct a poor man average trials
 200     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 201     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 202   }
 203   return kTRUE;
 204 }
 205
 206 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 207 {
 208
 209   //
 210   // Create the output container
 211   //
 212
 213
 214   // Connect the AOD
 215
 216
 217   MakeJetContainer();
 218
 219
 220   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 221
 222   OpenFile(1);
 223   if(!fHistList)fHistList = new TList();
 224
 225   Bool_t oldStatus = TH1::AddDirectoryStatus();
 226   TH1::AddDirectory(kFALSE);
 227
 228   //
 229   //  Histogram
 230
 231   const Int_t nBinPt = 100;
 232   Double_t binLimitsPt[nBinPt+1];
 233   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 234     if(iPt == 0){
 235       binLimitsPt[iPt] = 0.0;
 236     }
 237     else {// 1.0
 238       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 239     }
 240   }
 241
 242   const Int_t nBinPhi = 90;
 243   Double_t binLimitsPhi[nBinPhi+1];
 244   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 245     if(iPhi==0){
 246       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 247     }
 248     else{
 249       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 250     }
 251   }
 252
 253
 254
 255   const Int_t nBinEta = 40;
 256   Double_t binLimitsEta[nBinEta+1];
 257   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 258     if(iEta==0){
 259       binLimitsEta[iEta] = -2.0;
 260     }
 261     else{
 262       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 1/(Float_t)nBinEta + 0.1;
 263     }
 264   }
 265
 266   const Int_t nBinFrag = 25;
 267
 268
 269   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 270   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 271
 272   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 273   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 274
 275   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 276   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 277   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 278
 279   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 280   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 281
 282   fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 283   fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
 284   fh1SumPtTrackAreaRec = new TH1F("fh1SumPtTrackAreaRec","Sum Rec track P_T #eta <0.9 / 1.8 * 2 * 0.4*0.4;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
 285
 286   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 287   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 288   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 289   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 290   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 291
 292   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 293   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 294   //
 295
 296   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 297                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 298   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 299                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 300   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 301                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 302   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 303                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 304
 305   for(int ij = 0;ij < kMaxJets;++ij){
 306     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 307     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 308
 309     fh2PhiPt[ij] =  new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
 310                              nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 311
 312     fh2PhiEta[ij] =  new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
 313                               nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 314
 315
 316     fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
 317                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 318     fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 319                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 320
 321     fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
 322                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 323     fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 324                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 325   }
 326
 327
 328   const Int_t saveLevel = 3; // large save level more histos
 329   if(saveLevel>0){
 330     fHistList->Add(fh1Xsec);
 331     fHistList->Add(fh1Trials);
 332     fHistList->Add(fh1PtHard);
 333     fHistList->Add(fh1PtHardNoW);
 334     fHistList->Add(fh1PtHardTrials);
 335     if(fBranchGen.Length()>0){
 336       fHistList->Add(fh1NGenJets);
 337       fHistList->Add(fh1PtTracksGenIn);
 338     }
 339     fHistList->Add(fh1PtJetsRecIn);
 340     fHistList->Add(fh1PtJetsLeadingRecIn);
 341     fHistList->Add(fh1PtTracksRecIn);
 342     fHistList->Add(fh1PtTracksLeadingRecIn);
 343     fHistList->Add(fh1NRecJets);
 344     fHistList->Add(fh1PtTrackRec);
 345     fHistList->Add(fh1SumPtTrackRec);
 346     fHistList->Add(fh1SumPtTrackAreaRec);
 347     fHistList->Add(fh2NRecJetsPt);
 348     fHistList->Add(fh2NRecTracksPt);
 349     fHistList->Add(fh2JetsLeadingPhiEta );
 350     fHistList->Add(fh2JetsLeadingPhiPt );
 351     fHistList->Add(fh2TracksLeadingPhiEta);
 352     fHistList->Add(fh2TracksLeadingPhiPt);
 353     for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 354     for(int ij = 0;ij<kMaxJets;++ij){
 355       fHistList->Add( fh1PtRecIn[ij]);
 356       if(fBranchGen.Length()>0){
 357         fHistList->Add( fh1PtGenIn[ij]);
 358         fHistList->Add( fh2FragGen[ij]);
 359         fHistList->Add( fh2FragLnGen[ij]);
 360       }
 361       fHistList->Add( fh2PhiPt[ij]);
 362       fHistList->Add( fh2PhiEta[ij]);
 363       fHistList->Add( fh2FragRec[ij]);
 364       fHistList->Add( fh2FragLnRec[ij]);
 365     }
 366     fHistList->Add(fhnCorrelation);
 367   }
 368
 369   // =========== Switch on Sumw2 for all histos ===========
 370   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 371     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 372     if (h1){
 373       h1->Sumw2();
 374       continue;
 375     }
 376     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 377     if(hn)hn->Sumw2();
 378   }
 379   TH1::AddDirectory(oldStatus);
 380 }
 381
 382 void AliAnalysisTaskJetSpectrum2::Init()
 383 {
 384   //
 385   // Initialization
 386   //
 387
 388   Printf(">>> AnalysisTaskJetSpectrum2::Init() debug level %d\n",fDebug);
 389   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 390
 391 }
 392
 393 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 394 {
 395
 396   if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
 397     // no selection by the service task, we continue
 398     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 399     PostData(1, fHistList);
 400     return;
 401   }
 402
 403   //
 404   // Execute analysis for current event
 405   //
 406   AliESDEvent *fESD = 0;
 407   if(fUseAODJetInput){
 408     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 409     if(!fAOD){
 410       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
 411       return;
 412     }
 413     // fethc the header
 414   }
 415   else{
 416     //  assume that the AOD is in the general output...
 417     fAOD  = AODEvent();
 418     if(!fAOD){
 419       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 420       return;
 421     }
 422     if(fDebug>0){
 423       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 424     }
 425   }
 426
 427
 428
 429
 430   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 431
 432
 433   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 434
 435   if(!aodH){
 436     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 437     return;
 438   }
 439
 440   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 441   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 442   if(!aodRecJets){
 443     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 444     return;
 445   }
 446
 447   // ==== General variables needed
 448
 449
 450   // We use statice array, not to fragment the memory
 451   AliAODJet genJets[kMaxJets];
 452   Int_t nGenJets = 0;
 453   AliAODJet recJets[kMaxJets];
 454   Int_t nRecJets = 0;
 455   ///////////////////////////
 456
 457
 458   Double_t eventW = 1;
 459   Double_t ptHard = 0;
 460   Double_t nTrials = 1; // Trials for MC trigger
 461
 462   if(fUseExternalWeightOnly){
 463     eventW = fExternalWeight;
 464   }
 465
 466   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 467
 468   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 469   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 470     // this is the part we only use when we have MC information
 471     AliMCEvent* mcEvent = MCEvent();
 472     //    AliStack *pStack = 0;
 473     if(!mcEvent){
 474       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 475       return;
 476     }
 477     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 478     Int_t iCount = 0;
 479     if(pythiaGenHeader){
 480       nTrials = pythiaGenHeader->Trials();
 481       ptHard  = pythiaGenHeader->GetPtHard();
 482       int iProcessType = pythiaGenHeader->ProcessType();
 483       // 11 f+f -> f+f
 484       // 12 f+barf -> f+barf
 485       // 13 f+barf -> g+g
 486       // 28 f+g -> f+g
 487       // 53 g+g -> f+barf
 488       // 68 g+g -> g+g
 489       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 490       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 491
 492       // fetch the pythia generated jets only to be used here
 493       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 494       AliAODJet pythiaGenJets[kMaxJets];
 495       for(int ip = 0;ip < nPythiaGenJets;++ip){
 496         if(iCount>=kMaxJets)continue;
 497         Float_t p[4];
 498         pythiaGenHeader->TriggerJet(ip,p);
 499         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 500
 501         if(fLimitGenJetEta){
 502           if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 503              pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 504         }
 505
 506
 507         if(fBranchGen.Length()==0){
 508           // if we have MC particles and we do not read from the aod branch
 509           // use the pythia jets
 510           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 511         }
 512         iCount++;
 513       }
 514     }
 515     if(fBranchGen.Length()==0)nGenJets = iCount;
 516   }// (fAnalysisType&kMCESD)==kMCESD)
 517
 518
 519   // we simply fetch the tracks/mc particles as a list of AliVParticles
 520
 521   TList recParticles;
 522   TList genParticles;
 523
 524
 525
 526
 527   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 528   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 529   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 530   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 531
 532
 533   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 534   fh1PtHard->Fill(ptHard,eventW);
 535   fh1PtHardNoW->Fill(ptHard,1);
 536   fh1PtHardTrials->Fill(ptHard,nTrials);
 537
 538   // If we set a second branch for the input jets fetch this
 539   if(fBranchGen.Length()>0){
 540     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 541     if(aodGenJets){
 542       Int_t iCount = 0;
 543       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 544         if(iCount>=kMaxJets)continue;
 545         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 546         if(!tmp)continue;
 547         if(fLimitGenJetEta){
 548           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 549              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 550         }
 551         genJets[iCount] = *tmp;
 552         iCount++;
 553       }
 554       nGenJets = iCount;
 555     }
 556     else{
 557       if(fDebug>0)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 558       if(fDebug>1)fAOD->Print();
 559     }
 560   }
 561
 562   fh1NGenJets->Fill(nGenJets);
 563   // We do not want to exceed the maximum jet number
 564   nGenJets = TMath::Min(nGenJets,kMaxJets);
 565
 566   // Fetch the reconstructed jets...
 567
 568   nRecJets = aodRecJets->GetEntries();
 569
 570   nRecJets = aodRecJets->GetEntries();
 571   fh1NRecJets->Fill(nRecJets);
 572
 573   // Do something with the all rec jets
 574   Int_t nRecOver = nRecJets;
 575
 576   // check that the jets are sorted
 577   Float_t ptOld = 999999;
 578   for(int ir = 0;ir < nRecJets;ir++){
 579     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
 580     Float_t tmpPt = tmp->Pt();
 581     if(tmpPt>ptOld){
 582       Printf("%s:%d Jets Not Sorted!! %d:%.3E %d%.3E",(char*)__FILE__,__LINE__,ir,tmpPt,ir-1,ptOld);
 583     }
 584     ptOld = tmpPt;
 585   }
 586
 587
 588   if(nRecOver>0){
 589     TIterator *recIter = aodRecJets->MakeIterator();
 590     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
 591     Float_t pt = tmpRec->Pt();
 592     if(tmpRec){
 593       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 594         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 595         while(pt<ptCut&&tmpRec){
 596           nRecOver--;
 597           tmpRec = (AliAODJet*)(recIter->Next());
 598           if(tmpRec){
 599             pt = tmpRec->Pt();
 600           }
 601         }
 602         if(nRecOver<=0)break;
 603         fh2NRecJetsPt->Fill(ptCut,nRecOver);
 604       }
 605     }
 606     recIter->Reset();
 607     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
 608     Float_t phi = leading->Phi();
 609     if(phi<0)phi+=TMath::Pi()*2.;
 610     Float_t eta = leading->Eta();
 611      while((tmpRec = (AliAODJet*)(recIter->Next()))){
 612       Float_t tmpPt = tmpRec->Pt();
 613       fh1PtJetsRecIn->Fill(tmpPt);
 614       if(tmpRec==leading){
 615         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 616         continue;
 617       }
 618       // correlation
 619       Float_t tmpPhi =  tmpRec->Phi();
 620
 621       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 622       Float_t dPhi = phi - tmpRec->Phi();
 623       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 624       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 625       Float_t dEta = eta - tmpRec->Eta();
 626       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 627       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 628     }
 629     delete recIter;
 630   }
 631
 632   Int_t nTrackOver = recParticles.GetSize();
 633   // do the same for tracks and jets
 634   if(nTrackOver>0){
 635     TIterator *recIter = recParticles.MakeIterator();
 636     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 637     Float_t pt = tmpRec->Pt();
 638     //    Printf("Leading track p_t %3.3E",pt);
 639     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 640       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 641       while(pt<ptCut&&tmpRec){
 642         nTrackOver--;
 643         tmpRec = (AliVParticle*)(recIter->Next());
 644         if(tmpRec){
 645           pt = tmpRec->Pt();
 646         }
 647       }
 648       if(nTrackOver<=0)break;
 649       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 650     }
 651
 652     recIter->Reset();
 653     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 654     Float_t phi = leading->Phi();
 655     if(phi<0)phi+=TMath::Pi()*2.;
 656     Float_t eta = leading->Eta();
 657     while((tmpRec = (AliVParticle*)(recIter->Next()))){
 658       Float_t tmpPt = tmpRec->Pt();
 659       fh1PtTracksRecIn->Fill(tmpPt);
 660       if(tmpRec==leading){
 661         fh1PtTracksLeadingRecIn->Fill(tmpPt);
 662         continue;
 663       }
 664       // correlation
 665       Float_t tmpPhi =  tmpRec->Phi();
 666
 667       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 668       Float_t dPhi = phi - tmpRec->Phi();
 669       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 670       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 671       Float_t dEta = eta - tmpRec->Eta();
 672       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 673       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 674     }
 675     delete recIter;
 676   }
 677
 678   if(genParticles.GetSize()){
 679     TIterator *genIter = genParticles.MakeIterator();
 680     AliVParticle *tmpGen = 0;
 681     while((tmpGen = (AliVParticle*)(genIter->Next()))){
 682       if(TMath::Abs(tmpGen->Eta())<0.9){
 683         Float_t tmpPt = tmpGen->Pt();
 684         fh1PtTracksGenIn->Fill(tmpPt);
 685       }
 686     }
 687     delete genIter;
 688   }
 689
 690   fh1NRecJets->Fill(nRecJets);
 691   nRecJets = TMath::Min(nRecJets,kMaxJets);
 692
 693   for(int ir = 0;ir < nRecJets;++ir){
 694     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 695     if(!tmp)continue;
 696     recJets[ir] = *tmp;
 697   }
 698
 699
 700   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 701   // Relate the jets
 702   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 703   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 704
 705
 706   for(int i = 0;i<kMaxJets;++i){
 707     iGenIndex[i] = iRecIndex[i] = -1;
 708   }
 709
 710   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 711                                             iGenIndex,iRecIndex,fDebug);
 712   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 713
 714   if(fDebug){
 715     for(int i = 0;i<kMaxJets;++i){
 716       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 717       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 718     }
 719   }
 720
 721
 722
 723
 724   Double_t container[6];
 725
 726   // loop over generated jets
 727
 728   // radius; tmp, get from the jet header itself
 729   // at some point, how todeal woht FastJet on MC?
 730   Float_t radiusGen =0.4;
 731   Float_t radiusRec =0.4;
 732
 733   for(int ig = 0;ig < nGenJets;++ig){
 734     Double_t ptGen = genJets[ig].Pt();
 735     Double_t phiGen = genJets[ig].Phi();
 736     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 737     Double_t etaGen = genJets[ig].Eta();
 738
 739     container[3] = ptGen;
 740     container[4] = etaGen;
 741     container[5] = phiGen;
 742     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
 743     Int_t ir = iRecIndex[ig];
 744
 745     if(TMath::Abs(etaGen)<fRecEtaWindow){
 746       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
 747       fh1PtGenIn[ig]->Fill(ptGen,eventW);
 748       // fill the fragmentation function
 749       for(int it = 0;it<genParticles.GetEntries();++it){
 750         AliVParticle *part = (AliVParticle*)genParticles.At(it);
 751         if(genJets[ig].DeltaR(part)<radiusGen){
 752           Float_t z = part->Pt()/ptGen;
 753           Float_t lnz =  -1.*TMath::Log(z);
 754           fh2FragGen[ig]->Fill(z,ptGen,eventW);
 755           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
 756         }
 757       }
 758       if(ir>=0&&ir<nRecJets){
 759         fhnJetContainer[kStep3]->Fill(&container[3],eventW);
 760       }
 761     }
 762
 763     if(ir>=0&&ir<nRecJets){
 764       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
 765       Double_t etaRec = recJets[ir].Eta();
 766       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
 767     }
 768   }// loop over generated jets
 769
 770
 771   Float_t sumPt = 0;
 772   for(int it = 0;it<recParticles.GetEntries();++it){
 773     AliVParticle *part = (AliVParticle*)recParticles.At(it);
 774     // fill sum pt and P_t of all paritles
 775     if(TMath::Abs(part->Eta())<0.9){
 776       Float_t pt = part->Pt();
 777       fh1PtTrackRec->Fill(pt,eventW);
 778       sumPt += pt;
 779     }
 780   }
 781   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
 782   fh1SumPtTrackRec->Fill(sumPt,eventW);
 783
 784
 785   // loop over reconstructed jets
 786   for(int ir = 0;ir < nRecJets;++ir){
 787     Double_t etaRec = recJets[ir].Eta();
 788     Double_t ptRec = recJets[ir].Pt();
 789     Double_t phiRec = recJets[ir].Phi();
 790     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 791     container[0] = ptRec;
 792     container[1] = etaRec;
 793     container[2] = phiRec;
 794
 795     if(ptRec>10.&&fDebug>0){
 796       // need to cast to int, otherwise the printf overwrites
 797       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
 798       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
 799       fAOD->GetHeader()->Print();
 800       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
 801       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
 802         AliAODTrack *tr = fAOD->GetTrack(it);
 803         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 804         tr->Print();
 805         tr->Dump();
 806         if(fESD){
 807           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
 808           esdTr->Print("");
 809           esdTr->Dump();
 810         }
 811       }
 812     }
 813
 814
 815     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
 816     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 817     if(TMath::Abs(etaRec)<fRecEtaWindow){
 818       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
 819       fh1PtRecIn[ir]->Fill(ptRec,eventW);
 820       // fill the fragmentation function
 821       for(int it = 0;it<recParticles.GetEntries();++it){
 822         AliVParticle *part = (AliVParticle*)recParticles.At(it);
 823         Float_t eta = part->Eta();
 824         if(TMath::Abs(eta)<0.9){
 825           Float_t phi = part->Phi();
 826           if(phi<0)phi+=TMath::Pi()*2.;
 827           Float_t dPhi = phi - phiRec;
 828           Float_t dEta = eta - etaRec;
 829           if(dPhi<(-1.*TMath::Pi()))phiRec+=TMath::Pi()*2.;
 830           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
 831           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
 832         }
 833         if(recJets[ir].DeltaR(part)<radiusRec){
 834           Float_t z = part->Pt()/ptRec;
 835           Float_t lnz =  -1.*TMath::Log(z);
 836           fh2FragRec[ir]->Fill(z,ptRec,eventW);
 837           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
 838         }
 839       }
 840     }
 841
 842
 843     // Fill Correlation
 844     Int_t ig = iGenIndex[ir];
 845     if(ig>=0 && ig<nGenJets){
 846       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
 847       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 848       Double_t ptGen  = genJets[ig].Pt();
 849       Double_t phiGen = genJets[ig].Phi();
 850       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 851       Double_t etaGen = genJets[ig].Eta();
 852
 853       container[3] = ptGen;
 854       container[4] = etaGen;
 855       container[5] = phiGen;
 856
 857       //
 858       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
 859       //
 860
 861       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
 862       if(TMath::Abs(etaRec)<fRecEtaWindow){
 863         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
 864         fhnCorrelation->Fill(container);
 865       }// if etarec in window
 866
 867     }
 868     ////////////////////////////////////////////////////
 869     else{
 870
 871     }
 872   }// loop over reconstructed jets
 873
 874
 875   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 876   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 877   PostData(1, fHistList);
 878 }
 879
 880
 881 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
 882   //
 883   // Create the particle container for the correction framework manager and
 884   // link it
 885   //
 886   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 887   const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
 888   const Double_t kEtamin = -3.0, kEtamax = 3.0;
 889   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 890
 891   // can we neglect migration in eta and phi?
 892   // phi should be no problem since we cover full phi and are phi symmetric
 893   // eta migration is more difficult  due to needed acceptance correction
 894   // in limited eta range
 895
 896
 897   //arrays for the number of bins in each dimension
 898   Int_t iBin[kNvar];
 899   iBin[0] = 100; //bins in pt
 900   iBin[1] =  1; //bins in eta
 901   iBin[2] = 1; // bins in phi
 902
 903   //arrays for lower bounds :
 904   Double_t* binEdges[kNvar];
 905   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 906     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 907
 908   //values for bin lower bounds
 909   //  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
 910   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
 911   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 912   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 913
 914
 915   for(int i = 0;i < kMaxStep*2;++i){
 916     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
 917     for (int k=0; k<kNvar; k++) {
 918       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
 919     }
 920   }
 921   //create correlation matrix for unfolding
 922   Int_t thnDim[2*kNvar];
 923   for (int k=0; k<kNvar; k++) {
 924     //first half  : reconstructed
 925     //second half : MC
 926     thnDim[k]      = iBin[k];
 927     thnDim[k+kNvar] = iBin[k];
 928   }
 929
 930   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
 931   for (int k=0; k<kNvar; k++) {
 932     fhnCorrelation->SetBinEdges(k,binEdges[k]);
 933     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
 934   }
 935   fhnCorrelation->Sumw2();
 936
 937   // Add a histogram for Fake jets
 938   //  thnDim[3] = AliPID::kSPECIES;
 939   //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
 940   // for(Int_t idim = 0; idim < kNvar; idim++)
 941   //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
 942 }
 943
 944 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
 945 {
 946 // Terminate analysis
 947 //
 948     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
 949 }
 950
 951
 952 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
 953
 954   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 955
 956   Int_t iCount = 0;
 957   if(type==kTrackAOD){
 958     AliAODEvent *aod = 0;
 959     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 960     else aod = AODEvent();
 961     if(!aod){
 962       return iCount;
 963     }
 964     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 965       AliAODTrack *tr = aod->GetTrack(it);
 966       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 967       if(TMath::Abs(tr->Eta())>0.9)continue;
 968       list->Add(tr);
 969       iCount++;
 970     }
 971   }
 972   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 973     AliMCEvent* mcEvent = MCEvent();
 974     if(!mcEvent)return iCount;
 975     // we want to have alivpartilces so use get track
 976     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 977       if(!mcEvent->IsPhysicalPrimary(it))continue;
 978       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 979       if(type == kTrackKineAll){
 980         list->Add(part);
 981         iCount++;
 982       }
 983       else if(type == kTrackKineCharged){
 984         if(part->Particle()->GetPDG()->Charge()==0)continue;
 985         list->Add(part);
 986         iCount++;
 987       }
 988     }
 989   }
 990   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 991     AliAODEvent *aod = 0;
 992     if(fUseAODMCInput)dynamic_cast<AliAODEvent*>(InputEvent());
 993     else aod = AODEvent();
 994     if(!aod)return iCount;
 995     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 996     if(!tca)return iCount;
 997     for(int it = 0;it < tca->GetEntriesFast();++it){
 998       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 999       if(!part->IsPhysicalPrimary())continue;
1000       if(type == kTrackAODMCAll){
1001         list->Add(part);
1002         iCount++;
1003       }
1004       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1005         if(part->Charge()==0)continue;
1006         if(kTrackAODMCCharged){
1007           list->Add(part);
1008         }
1009         else {
1010           if(TMath::Abs(part->Eta())>0.9)continue;
1011           list->Add(part);
1012         }
1013         iCount++;
1014       }
1015     }
1016   }// AODMCparticle
1017   list->Sort();
1018   return iCount;
1019
1020 }