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   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 391
 392 }
 393
 394 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 395 {
 396
 397   if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
 398     // no selection by the service task, we continue
 399     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 400     PostData(1, fHistList);
 401     return;
 402   }
 403
 404   //
 405   // Execute analysis for current event
 406   //
 407   AliESDEvent *fESD = 0;
 408   if(fUseAODJetInput){
 409     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 410     if(!fAOD){
 411       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
 412       return;
 413     }
 414     // fethc the header
 415   }
 416   else{
 417     //  assume that the AOD is in the general output...
 418     fAOD  = AODEvent();
 419     if(!fAOD){
 420       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 421       return;
 422     }
 423     if(fDebug>0){
 424       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 425     }
 426   }
 427
 428
 429
 430
 431   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 432
 433
 434   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 435
 436   if(!aodH){
 437     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 438     return;
 439   }
 440
 441   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 442   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 443   if(!aodRecJets){
 444     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 445     return;
 446   }
 447
 448   // ==== General variables needed
 449
 450
 451   // We use statice array, not to fragment the memory
 452   AliAODJet genJets[kMaxJets];
 453   Int_t nGenJets = 0;
 454   AliAODJet recJets[kMaxJets];
 455   Int_t nRecJets = 0;
 456   ///////////////////////////
 457
 458
 459   Double_t eventW = 1;
 460   Double_t ptHard = 0;
 461   Double_t nTrials = 1; // Trials for MC trigger
 462
 463   if(fUseExternalWeightOnly){
 464     eventW = fExternalWeight;
 465   }
 466
 467   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 468   //  if(fDebug>0)aodH->SetFillAOD(kFALSE);
 469   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 470   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 471     // this is the part we only use when we have MC information
 472     AliMCEvent* mcEvent = MCEvent();
 473     //    AliStack *pStack = 0;
 474     if(!mcEvent){
 475       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 476       return;
 477     }
 478     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 479     Int_t iCount = 0;
 480     if(pythiaGenHeader){
 481       nTrials = pythiaGenHeader->Trials();
 482       ptHard  = pythiaGenHeader->GetPtHard();
 483       int iProcessType = pythiaGenHeader->ProcessType();
 484       // 11 f+f -> f+f
 485       // 12 f+barf -> f+barf
 486       // 13 f+barf -> g+g
 487       // 28 f+g -> f+g
 488       // 53 g+g -> f+barf
 489       // 68 g+g -> g+g
 490       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 491       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 492
 493       // fetch the pythia generated jets only to be used here
 494       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 495       AliAODJet pythiaGenJets[kMaxJets];
 496       for(int ip = 0;ip < nPythiaGenJets;++ip){
 497         if(iCount>=kMaxJets)continue;
 498         Float_t p[4];
 499         pythiaGenHeader->TriggerJet(ip,p);
 500         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 501
 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         /*
 550         if(fLimitGenJetEta){
 551           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 552              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 553         }
 554         */
 555         genJets[iCount] = *tmp;
 556         iCount++;
 557       }
 558       nGenJets = iCount;
 559     }
 560     else{
 561       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 562       if(fDebug>2)fAOD->Print();
 563     }
 564   }
 565
 566   fh1NGenJets->Fill(nGenJets);
 567   // We do not want to exceed the maximum jet number
 568   nGenJets = TMath::Min(nGenJets,kMaxJets);
 569
 570   // Fetch the reconstructed jets...
 571
 572   nRecJets = aodRecJets->GetEntries();
 573
 574   nRecJets = aodRecJets->GetEntries();
 575   fh1NRecJets->Fill(nRecJets);
 576
 577   // Do something with the all rec jets
 578   Int_t nRecOver = nRecJets;
 579
 580   // check that the jets are sorted
 581   Float_t ptOld = 999999;
 582   for(int ir = 0;ir < nRecJets;ir++){
 583     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
 584     Float_t tmpPt = tmp->Pt();
 585     if(tmpPt>ptOld){
 586       Printf("%s:%d Jets Not Sorted!! %d:%.3E %d%.3E",(char*)__FILE__,__LINE__,ir,tmpPt,ir-1,ptOld);
 587     }
 588     ptOld = tmpPt;
 589   }
 590
 591
 592   if(nRecOver>0){
 593     TIterator *recIter = aodRecJets->MakeIterator();
 594     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
 595     Float_t pt = tmpRec->Pt();
 596     if(tmpRec){
 597       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 598         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 599         while(pt<ptCut&&tmpRec){
 600           nRecOver--;
 601           tmpRec = (AliAODJet*)(recIter->Next());
 602           if(tmpRec){
 603             pt = tmpRec->Pt();
 604           }
 605         }
 606         if(nRecOver<=0)break;
 607         fh2NRecJetsPt->Fill(ptCut,nRecOver);
 608       }
 609     }
 610     recIter->Reset();
 611     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
 612     Float_t phi = leading->Phi();
 613     if(phi<0)phi+=TMath::Pi()*2.;
 614     Float_t eta = leading->Eta();
 615     pt = leading->Pt();
 616      while((tmpRec = (AliAODJet*)(recIter->Next()))){
 617       Float_t tmpPt = tmpRec->Pt();
 618       fh1PtJetsRecIn->Fill(tmpPt);
 619       if(tmpRec==leading){
 620         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 621         continue;
 622       }
 623       // correlation
 624       Float_t tmpPhi =  tmpRec->Phi();
 625
 626       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 627       Float_t dPhi = phi - tmpRec->Phi();
 628       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 629       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 630       Float_t dEta = eta - tmpRec->Eta();
 631       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 632       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 633     }
 634     delete recIter;
 635   }
 636
 637   Int_t nTrackOver = recParticles.GetSize();
 638   // do the same for tracks and jets
 639   if(nTrackOver>0){
 640     TIterator *recIter = recParticles.MakeIterator();
 641     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 642     Float_t pt = tmpRec->Pt();
 643     //    Printf("Leading track p_t %3.3E",pt);
 644     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 645       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 646       while(pt<ptCut&&tmpRec){
 647         nTrackOver--;
 648         tmpRec = (AliVParticle*)(recIter->Next());
 649         if(tmpRec){
 650           pt = tmpRec->Pt();
 651         }
 652       }
 653       if(nTrackOver<=0)break;
 654       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 655     }
 656
 657     recIter->Reset();
 658     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 659     Float_t phi = leading->Phi();
 660     if(phi<0)phi+=TMath::Pi()*2.;
 661     Float_t eta = leading->Eta();
 662     pt  = leading->Pt();
 663     while((tmpRec = (AliVParticle*)(recIter->Next()))){
 664       Float_t tmpPt = tmpRec->Pt();
 665       fh1PtTracksRecIn->Fill(tmpPt);
 666       if(tmpRec==leading){
 667         fh1PtTracksLeadingRecIn->Fill(tmpPt);
 668         continue;
 669       }
 670       // correlation
 671       Float_t tmpPhi =  tmpRec->Phi();
 672
 673       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 674       Float_t dPhi = phi - tmpRec->Phi();
 675       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 676       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 677       Float_t dEta = eta - tmpRec->Eta();
 678       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 679       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 680     }
 681     delete recIter;
 682   }
 683
 684   if(genParticles.GetSize()){
 685     TIterator *genIter = genParticles.MakeIterator();
 686     AliVParticle *tmpGen = 0;
 687     while((tmpGen = (AliVParticle*)(genIter->Next()))){
 688       if(TMath::Abs(tmpGen->Eta())<0.9){
 689         Float_t tmpPt = tmpGen->Pt();
 690         fh1PtTracksGenIn->Fill(tmpPt);
 691       }
 692     }
 693     delete genIter;
 694   }
 695
 696   fh1NRecJets->Fill(nRecJets);
 697   nRecJets = TMath::Min(nRecJets,kMaxJets);
 698
 699   for(int ir = 0;ir < nRecJets;++ir){
 700     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 701     if(!tmp)continue;
 702     recJets[ir] = *tmp;
 703   }
 704
 705
 706   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 707   // Relate the jets
 708   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 709   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 710
 711
 712   for(int i = 0;i<kMaxJets;++i){
 713     iGenIndex[i] = iRecIndex[i] = -1;
 714   }
 715
 716   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 717                                             iGenIndex,iRecIndex,fDebug);
 718   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 719
 720   if(fDebug){
 721     for(int i = 0;i<kMaxJets;++i){
 722       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 723       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 724     }
 725   }
 726
 727
 728
 729
 730   Double_t container[6];
 731
 732   // loop over generated jets
 733
 734   // radius; tmp, get from the jet header itself
 735   // at some point, how todeal woht FastJet on MC?
 736   Float_t radiusGen =0.4;
 737   Float_t radiusRec =0.4;
 738
 739   for(int ig = 0;ig < nGenJets;++ig){
 740     Double_t ptGen = genJets[ig].Pt();
 741     Double_t phiGen = genJets[ig].Phi();
 742     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 743     Double_t etaGen = genJets[ig].Eta();
 744
 745     container[3] = ptGen;
 746     container[4] = etaGen;
 747     container[5] = phiGen;
 748     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
 749     Int_t ir = iRecIndex[ig];
 750
 751     if(TMath::Abs(etaGen)<fRecEtaWindow){
 752       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
 753       fh1PtGenIn[ig]->Fill(ptGen,eventW);
 754       // fill the fragmentation function
 755       for(int it = 0;it<genParticles.GetEntries();++it){
 756         AliVParticle *part = (AliVParticle*)genParticles.At(it);
 757         if(genJets[ig].DeltaR(part)<radiusGen){
 758           Float_t z = part->Pt()/ptGen;
 759           Float_t lnz =  -1.*TMath::Log(z);
 760           fh2FragGen[ig]->Fill(z,ptGen,eventW);
 761           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
 762         }
 763       }
 764       if(ir>=0&&ir<nRecJets){
 765         fhnJetContainer[kStep3]->Fill(&container[3],eventW);
 766       }
 767     }
 768
 769     if(ir>=0&&ir<nRecJets){
 770       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
 771       Double_t etaRec = recJets[ir].Eta();
 772       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
 773     }
 774   }// loop over generated jets
 775
 776
 777   Float_t sumPt = 0;
 778   for(int it = 0;it<recParticles.GetEntries();++it){
 779     AliVParticle *part = (AliVParticle*)recParticles.At(it);
 780     // fill sum pt and P_t of all paritles
 781     if(TMath::Abs(part->Eta())<0.9){
 782       Float_t pt = part->Pt();
 783       fh1PtTrackRec->Fill(pt,eventW);
 784       sumPt += pt;
 785     }
 786   }
 787   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
 788   fh1SumPtTrackRec->Fill(sumPt,eventW);
 789
 790
 791   // loop over reconstructed jets
 792   for(int ir = 0;ir < nRecJets;++ir){
 793     Double_t etaRec = recJets[ir].Eta();
 794     Double_t ptRec = recJets[ir].Pt();
 795     Double_t phiRec = recJets[ir].Phi();
 796     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 797     container[0] = ptRec;
 798     container[1] = etaRec;
 799     container[2] = phiRec;
 800
 801     if(ptRec>30.&&fDebug>0){
 802       // need to cast to int, otherwise the printf overwrites
 803       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
 804       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
 805       if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
 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 }