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