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   fHistList->SetOwner(kTRUE);
 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         /*
 504         if(fLimitGenJetEta){
 505           if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 506              pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 507         }
 508         */
 509
 510         if(fBranchGen.Length()==0){
 511           // if we have MC particles and we do not read from the aod branch
 512           // use the pythia jets
 513           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 514         }
 515         iCount++;
 516       }
 517     }
 518     if(fBranchGen.Length()==0)nGenJets = iCount;
 519   }// (fAnalysisType&kMCESD)==kMCESD)
 520
 521
 522   // we simply fetch the tracks/mc particles as a list of AliVParticles
 523
 524   TList recParticles;
 525   TList genParticles;
 526
 527
 528
 529
 530   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 531   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 532   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 533   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 534
 535
 536   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 537   fh1PtHard->Fill(ptHard,eventW);
 538   fh1PtHardNoW->Fill(ptHard,1);
 539   fh1PtHardTrials->Fill(ptHard,nTrials);
 540
 541   // If we set a second branch for the input jets fetch this
 542   if(fBranchGen.Length()>0){
 543     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 544     if(aodGenJets){
 545       Int_t iCount = 0;
 546       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 547         if(iCount>=kMaxJets)continue;
 548         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 549         if(!tmp)continue;
 550         /*
 551         if(fLimitGenJetEta){
 552           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 553              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 554         }
 555         */
 556         genJets[iCount] = *tmp;
 557         iCount++;
 558       }
 559       nGenJets = iCount;
 560     }
 561     else{
 562       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 563       if(fDebug>2)fAOD->Print();
 564     }
 565   }
 566
 567   fh1NGenJets->Fill(nGenJets);
 568   // We do not want to exceed the maximum jet number
 569   nGenJets = TMath::Min(nGenJets,kMaxJets);
 570
 571   // Fetch the reconstructed jets...
 572
 573   nRecJets = aodRecJets->GetEntries();
 574
 575   nRecJets = aodRecJets->GetEntries();
 576   fh1NRecJets->Fill(nRecJets);
 577
 578   // Do something with the all rec jets
 579   Int_t nRecOver = nRecJets;
 580
 581   // check that the jets are sorted
 582   Float_t ptOld = 999999;
 583   for(int ir = 0;ir < nRecJets;ir++){
 584     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
 585     Float_t tmpPt = tmp->Pt();
 586     if(tmpPt>ptOld){
 587       Printf("%s:%d Jets Not Sorted!! %d:%.3E %d%.3E",(char*)__FILE__,__LINE__,ir,tmpPt,ir-1,ptOld);
 588     }
 589     ptOld = tmpPt;
 590   }
 591
 592
 593   if(nRecOver>0){
 594     TIterator *recIter = aodRecJets->MakeIterator();
 595     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
 596     Float_t pt = tmpRec->Pt();
 597     if(tmpRec){
 598       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 599         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 600         while(pt<ptCut&&tmpRec){
 601           nRecOver--;
 602           tmpRec = (AliAODJet*)(recIter->Next());
 603           if(tmpRec){
 604             pt = tmpRec->Pt();
 605           }
 606         }
 607         if(nRecOver<=0)break;
 608         fh2NRecJetsPt->Fill(ptCut,nRecOver);
 609       }
 610     }
 611     recIter->Reset();
 612     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
 613     Float_t phi = leading->Phi();
 614     if(phi<0)phi+=TMath::Pi()*2.;
 615     Float_t eta = leading->Eta();
 616     pt = leading->Pt();
 617      while((tmpRec = (AliAODJet*)(recIter->Next()))){
 618       Float_t tmpPt = tmpRec->Pt();
 619       fh1PtJetsRecIn->Fill(tmpPt);
 620       if(tmpRec==leading){
 621         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 622         continue;
 623       }
 624       // correlation
 625       Float_t tmpPhi =  tmpRec->Phi();
 626
 627       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 628       Float_t dPhi = phi - tmpRec->Phi();
 629       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 630       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 631       Float_t dEta = eta - tmpRec->Eta();
 632       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 633       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 634     }
 635     delete recIter;
 636   }
 637
 638   Int_t nTrackOver = recParticles.GetSize();
 639   // do the same for tracks and jets
 640   if(nTrackOver>0){
 641     TIterator *recIter = recParticles.MakeIterator();
 642     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 643     Float_t pt = tmpRec->Pt();
 644     //    Printf("Leading track p_t %3.3E",pt);
 645     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 646       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 647       while(pt<ptCut&&tmpRec){
 648         nTrackOver--;
 649         tmpRec = (AliVParticle*)(recIter->Next());
 650         if(tmpRec){
 651           pt = tmpRec->Pt();
 652         }
 653       }
 654       if(nTrackOver<=0)break;
 655       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 656     }
 657
 658     recIter->Reset();
 659     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 660     Float_t phi = leading->Phi();
 661     if(phi<0)phi+=TMath::Pi()*2.;
 662     Float_t eta = leading->Eta();
 663     pt  = leading->Pt();
 664     while((tmpRec = (AliVParticle*)(recIter->Next()))){
 665       Float_t tmpPt = tmpRec->Pt();
 666       fh1PtTracksRecIn->Fill(tmpPt);
 667       if(tmpRec==leading){
 668         fh1PtTracksLeadingRecIn->Fill(tmpPt);
 669         continue;
 670       }
 671       // correlation
 672       Float_t tmpPhi =  tmpRec->Phi();
 673
 674       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 675       Float_t dPhi = phi - tmpRec->Phi();
 676       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 677       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 678       Float_t dEta = eta - tmpRec->Eta();
 679       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 680       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 681     }
 682     delete recIter;
 683   }
 684
 685   if(genParticles.GetSize()){
 686     TIterator *genIter = genParticles.MakeIterator();
 687     AliVParticle *tmpGen = 0;
 688     while((tmpGen = (AliVParticle*)(genIter->Next()))){
 689       if(TMath::Abs(tmpGen->Eta())<0.9){
 690         Float_t tmpPt = tmpGen->Pt();
 691         fh1PtTracksGenIn->Fill(tmpPt);
 692       }
 693     }
 694     delete genIter;
 695   }
 696
 697   fh1NRecJets->Fill(nRecJets);
 698   nRecJets = TMath::Min(nRecJets,kMaxJets);
 699
 700   for(int ir = 0;ir < nRecJets;++ir){
 701     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 702     if(!tmp)continue;
 703     recJets[ir] = *tmp;
 704   }
 705
 706
 707   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 708   // Relate the jets
 709   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 710   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 711
 712
 713   for(int i = 0;i<kMaxJets;++i){
 714     iGenIndex[i] = iRecIndex[i] = -1;
 715   }
 716
 717   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 718                                             iGenIndex,iRecIndex,fDebug);
 719   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 720
 721   if(fDebug){
 722     for(int i = 0;i<kMaxJets;++i){
 723       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 724       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 725     }
 726   }
 727
 728
 729
 730
 731   Double_t container[6];
 732
 733   // loop over generated jets
 734
 735   // radius; tmp, get from the jet header itself
 736   // at some point, how todeal woht FastJet on MC?
 737   Float_t radiusGen =0.4;
 738   Float_t radiusRec =0.4;
 739
 740   for(int ig = 0;ig < nGenJets;++ig){
 741     Double_t ptGen = genJets[ig].Pt();
 742     Double_t phiGen = genJets[ig].Phi();
 743     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 744     Double_t etaGen = genJets[ig].Eta();
 745
 746     container[3] = ptGen;
 747     container[4] = etaGen;
 748     container[5] = phiGen;
 749     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
 750     Int_t ir = iRecIndex[ig];
 751
 752     if(TMath::Abs(etaGen)<fRecEtaWindow){
 753       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
 754       fh1PtGenIn[ig]->Fill(ptGen,eventW);
 755       // fill the fragmentation function
 756       for(int it = 0;it<genParticles.GetEntries();++it){
 757         AliVParticle *part = (AliVParticle*)genParticles.At(it);
 758         if(genJets[ig].DeltaR(part)<radiusGen){
 759           Float_t z = part->Pt()/ptGen;
 760           Float_t lnz =  -1.*TMath::Log(z);
 761           fh2FragGen[ig]->Fill(z,ptGen,eventW);
 762           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
 763         }
 764       }
 765       if(ir>=0&&ir<nRecJets){
 766         fhnJetContainer[kStep3]->Fill(&container[3],eventW);
 767       }
 768     }
 769
 770     if(ir>=0&&ir<nRecJets){
 771       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
 772       Double_t etaRec = recJets[ir].Eta();
 773       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
 774     }
 775   }// loop over generated jets
 776
 777
 778   Float_t sumPt = 0;
 779   for(int it = 0;it<recParticles.GetEntries();++it){
 780     AliVParticle *part = (AliVParticle*)recParticles.At(it);
 781     // fill sum pt and P_t of all paritles
 782     if(TMath::Abs(part->Eta())<0.9){
 783       Float_t pt = part->Pt();
 784       fh1PtTrackRec->Fill(pt,eventW);
 785       sumPt += pt;
 786     }
 787   }
 788   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
 789   fh1SumPtTrackRec->Fill(sumPt,eventW);
 790
 791
 792   // loop over reconstructed jets
 793   for(int ir = 0;ir < nRecJets;++ir){
 794     Double_t etaRec = recJets[ir].Eta();
 795     Double_t ptRec = recJets[ir].Pt();
 796     Double_t phiRec = recJets[ir].Phi();
 797     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 798     container[0] = ptRec;
 799     container[1] = etaRec;
 800     container[2] = phiRec;
 801
 802     if(ptRec>10.&&fDebug>0){
 803       // need to cast to int, otherwise the printf overwrites
 804       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
 805       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
 806       //  aodH->SetFillAOD(kTRUE);
 807       fAOD->GetHeader()->Print();
 808       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
 809       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
 810         AliAODTrack *tr = fAOD->GetTrack(it);
 811         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 812         tr->Print();
 813         //      tr->Dump();
 814         if(fESD){
 815           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
 816           esdTr->Print("");
 817           esdTr->Dump();
 818         }
 819       }
 820     }
 821
 822
 823     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
 824     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 825     if(TMath::Abs(etaRec)<fRecEtaWindow){
 826       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
 827       fh1PtRecIn[ir]->Fill(ptRec,eventW);
 828       // fill the fragmentation function
 829       for(int it = 0;it<recParticles.GetEntries();++it){
 830         AliVParticle *part = (AliVParticle*)recParticles.At(it);
 831         Float_t eta = part->Eta();
 832         if(TMath::Abs(eta)<0.9){
 833           Float_t phi = part->Phi();
 834           if(phi<0)phi+=TMath::Pi()*2.;
 835           Float_t dPhi = phi - phiRec;
 836           Float_t dEta = eta - etaRec;
 837           if(dPhi<(-1.*TMath::Pi()))phiRec+=TMath::Pi()*2.;
 838           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
 839           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
 840         }
 841         if(recJets[ir].DeltaR(part)<radiusRec){
 842           Float_t z = part->Pt()/ptRec;
 843           Float_t lnz =  -1.*TMath::Log(z);
 844           fh2FragRec[ir]->Fill(z,ptRec,eventW);
 845           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
 846         }
 847       }
 848     }
 849
 850
 851     // Fill Correlation
 852     Int_t ig = iGenIndex[ir];
 853     if(ig>=0 && ig<nGenJets){
 854       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
 855       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 856       Double_t ptGen  = genJets[ig].Pt();
 857       Double_t phiGen = genJets[ig].Phi();
 858       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 859       Double_t etaGen = genJets[ig].Eta();
 860
 861       container[3] = ptGen;
 862       container[4] = etaGen;
 863       container[5] = phiGen;
 864
 865       //
 866       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
 867       //
 868
 869       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
 870       if(TMath::Abs(etaRec)<fRecEtaWindow){
 871         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
 872         fhnCorrelation->Fill(container);
 873       }// if etarec in window
 874
 875     }
 876     ////////////////////////////////////////////////////
 877     else{
 878
 879     }
 880   }// loop over reconstructed jets
 881
 882
 883   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 884   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 885   PostData(1, fHistList);
 886 }
 887
 888
 889 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
 890   //
 891   // Create the particle container for the correction framework manager and
 892   // link it
 893   //
 894   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 895   const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
 896   const Double_t kEtamin = -3.0, kEtamax = 3.0;
 897   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 898
 899   // can we neglect migration in eta and phi?
 900   // phi should be no problem since we cover full phi and are phi symmetric
 901   // eta migration is more difficult  due to needed acceptance correction
 902   // in limited eta range
 903
 904
 905   //arrays for the number of bins in each dimension
 906   Int_t iBin[kNvar];
 907   iBin[0] = 100; //bins in pt
 908   iBin[1] =  1; //bins in eta
 909   iBin[2] = 1; // bins in phi
 910
 911   //arrays for lower bounds :
 912   Double_t* binEdges[kNvar];
 913   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 914     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 915
 916   //values for bin lower bounds
 917   //  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);
 918   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
 919   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 920   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 921
 922
 923   for(int i = 0;i < kMaxStep*2;++i){
 924     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
 925     for (int k=0; k<kNvar; k++) {
 926       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
 927     }
 928   }
 929   //create correlation matrix for unfolding
 930   Int_t thnDim[2*kNvar];
 931   for (int k=0; k<kNvar; k++) {
 932     //first half  : reconstructed
 933     //second half : MC
 934     thnDim[k]      = iBin[k];
 935     thnDim[k+kNvar] = iBin[k];
 936   }
 937
 938   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
 939   for (int k=0; k<kNvar; k++) {
 940     fhnCorrelation->SetBinEdges(k,binEdges[k]);
 941     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
 942   }
 943   fhnCorrelation->Sumw2();
 944
 945   // Add a histogram for Fake jets
 946   //  thnDim[3] = AliPID::kSPECIES;
 947   //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
 948   // for(Int_t idim = 0; idim < kNvar; idim++)
 949   //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
 950   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 951     delete [] binEdges[ivar];
 952
 953 }
 954
 955 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
 956 {
 957 // Terminate analysis
 958 //
 959     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
 960 }
 961
 962
 963 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
 964
 965   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 966
 967   Int_t iCount = 0;
 968   if(type==kTrackAOD){
 969     AliAODEvent *aod = 0;
 970     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 971     else aod = AODEvent();
 972     if(!aod){
 973       return iCount;
 974     }
 975     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 976       AliAODTrack *tr = aod->GetTrack(it);
 977       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 978       if(TMath::Abs(tr->Eta())>0.9)continue;
 979       if(fDebug>0){
 980         if(tr->Pt()>20){
 981           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
 982           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
 983           tr->Print();
 984           //    tr->Dump();
 985           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 986           if(fESD){
 987             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
 988             esdTr->Print("");
 989             esdTr->Dump();
 990           }
 991         }
 992       }
 993       list->Add(tr);
 994       iCount++;
 995     }
 996   }
 997   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 998     AliMCEvent* mcEvent = MCEvent();
 999     if(!mcEvent)return iCount;
1000     // we want to have alivpartilces so use get track
1001     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1002       if(!mcEvent->IsPhysicalPrimary(it))continue;
1003       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1004       if(type == kTrackKineAll){
1005         list->Add(part);
1006         iCount++;
1007       }
1008       else if(type == kTrackKineCharged){
1009         if(part->Particle()->GetPDG()->Charge()==0)continue;
1010         list->Add(part);
1011         iCount++;
1012       }
1013     }
1014   }
1015   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1016     AliAODEvent *aod = 0;
1017     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1018     else aod = AODEvent();
1019     if(!aod)return iCount;
1020     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1021     if(!tca)return iCount;
1022     for(int it = 0;it < tca->GetEntriesFast();++it){
1023       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1024       if(!part->IsPhysicalPrimary())continue;
1025       if(type == kTrackAODMCAll){
1026         list->Add(part);
1027         iCount++;
1028       }
1029       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1030         if(part->Charge()==0)continue;
1031         if(kTrackAODMCCharged){
1032           list->Add(part);
1033         }
1034         else {
1035           if(TMath::Abs(part->Eta())>0.9)continue;
1036           list->Add(part);
1037         }
1038         iCount++;
1039       }
1040     }
1041   }// AODMCparticle
1042   list->Sort();
1043   return iCount;
1044
1045 }