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>10.&&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(),fInputHandler->GetTree()->GetReadEntry());
 805       //  aodH->SetFillAOD(kTRUE);
 806       fAOD->GetHeader()->Print();
 807       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
 808       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
 809         AliAODTrack *tr = fAOD->GetTrack(it);
 810         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 811         tr->Print();
 812         //      tr->Dump();
 813         if(fESD){
 814           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
 815           esdTr->Print("");
 816           // esdTr->Dump();
 817         }
 818       }
 819     }
 820
 821
 822     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
 823     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 824     if(TMath::Abs(etaRec)<fRecEtaWindow){
 825       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
 826       fh1PtRecIn[ir]->Fill(ptRec,eventW);
 827       // fill the fragmentation function
 828       for(int it = 0;it<recParticles.GetEntries();++it){
 829         AliVParticle *part = (AliVParticle*)recParticles.At(it);
 830         Float_t eta = part->Eta();
 831         if(TMath::Abs(eta)<0.9){
 832           Float_t phi = part->Phi();
 833           if(phi<0)phi+=TMath::Pi()*2.;
 834           Float_t dPhi = phi - phiRec;
 835           Float_t dEta = eta - etaRec;
 836           if(dPhi<(-1.*TMath::Pi()))phiRec+=TMath::Pi()*2.;
 837           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
 838           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
 839         }
 840         if(recJets[ir].DeltaR(part)<radiusRec){
 841           Float_t z = part->Pt()/ptRec;
 842           Float_t lnz =  -1.*TMath::Log(z);
 843           fh2FragRec[ir]->Fill(z,ptRec,eventW);
 844           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
 845         }
 846       }
 847     }
 848
 849
 850     // Fill Correlation
 851     Int_t ig = iGenIndex[ir];
 852     if(ig>=0 && ig<nGenJets){
 853       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
 854       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 855       Double_t ptGen  = genJets[ig].Pt();
 856       Double_t phiGen = genJets[ig].Phi();
 857       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 858       Double_t etaGen = genJets[ig].Eta();
 859
 860       container[3] = ptGen;
 861       container[4] = etaGen;
 862       container[5] = phiGen;
 863
 864       //
 865       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
 866       //
 867
 868       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
 869       if(TMath::Abs(etaRec)<fRecEtaWindow){
 870         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
 871         fhnCorrelation->Fill(container);
 872       }// if etarec in window
 873
 874     }
 875     ////////////////////////////////////////////////////
 876     else{
 877
 878     }
 879   }// loop over reconstructed jets
 880
 881
 882   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 883   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 884   PostData(1, fHistList);
 885 }
 886
 887
 888 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
 889   //
 890   // Create the particle container for the correction framework manager and
 891   // link it
 892   //
 893   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 894   const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
 895   const Double_t kEtamin = -3.0, kEtamax = 3.0;
 896   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 897
 898   // can we neglect migration in eta and phi?
 899   // phi should be no problem since we cover full phi and are phi symmetric
 900   // eta migration is more difficult  due to needed acceptance correction
 901   // in limited eta range
 902
 903
 904   //arrays for the number of bins in each dimension
 905   Int_t iBin[kNvar];
 906   iBin[0] = 100; //bins in pt
 907   iBin[1] =  1; //bins in eta
 908   iBin[2] = 1; // bins in phi
 909
 910   //arrays for lower bounds :
 911   Double_t* binEdges[kNvar];
 912   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 913     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 914
 915   //values for bin lower bounds
 916   //  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);
 917   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
 918   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 919   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 920
 921
 922   for(int i = 0;i < kMaxStep*2;++i){
 923     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
 924     for (int k=0; k<kNvar; k++) {
 925       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
 926     }
 927   }
 928   //create correlation matrix for unfolding
 929   Int_t thnDim[2*kNvar];
 930   for (int k=0; k<kNvar; k++) {
 931     //first half  : reconstructed
 932     //second half : MC
 933     thnDim[k]      = iBin[k];
 934     thnDim[k+kNvar] = iBin[k];
 935   }
 936
 937   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
 938   for (int k=0; k<kNvar; k++) {
 939     fhnCorrelation->SetBinEdges(k,binEdges[k]);
 940     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
 941   }
 942   fhnCorrelation->Sumw2();
 943
 944   // Add a histogram for Fake jets
 945   //  thnDim[3] = AliPID::kSPECIES;
 946   //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
 947   // for(Int_t idim = 0; idim < kNvar; idim++)
 948   //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
 949   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 950     delete [] binEdges[ivar];
 951
 952 }
 953
 954 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
 955 {
 956 // Terminate analysis
 957 //
 958     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
 959 }
 960
 961
 962 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
 963
 964   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 965
 966   Int_t iCount = 0;
 967   if(type==kTrackAOD){
 968     AliAODEvent *aod = 0;
 969     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 970     else aod = AODEvent();
 971     if(!aod){
 972       return iCount;
 973     }
 974     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 975       AliAODTrack *tr = aod->GetTrack(it);
 976       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 977       if(TMath::Abs(tr->Eta())>0.9)continue;
 978       if(fDebug>0){
 979         if(tr->Pt()>20){
 980           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
 981           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
 982           tr->Print();
 983           //    tr->Dump();
 984           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 985           if(fESD){
 986             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
 987             esdTr->Print("");
 988             // esdTr->Dump();
 989           }
 990         }
 991       }
 992       list->Add(tr);
 993       iCount++;
 994     }
 995   }
 996   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 997     AliMCEvent* mcEvent = MCEvent();
 998     if(!mcEvent)return iCount;
 999     // we want to have alivpartilces so use get track
1000     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1001       if(!mcEvent->IsPhysicalPrimary(it))continue;
1002       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1003       if(type == kTrackKineAll){
1004         list->Add(part);
1005         iCount++;
1006       }
1007       else if(type == kTrackKineCharged){
1008         if(part->Particle()->GetPDG()->Charge()==0)continue;
1009         list->Add(part);
1010         iCount++;
1011       }
1012     }
1013   }
1014   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1015     AliAODEvent *aod = 0;
1016     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1017     else aod = AODEvent();
1018     if(!aod)return iCount;
1019     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1020     if(!tca)return iCount;
1021     for(int it = 0;it < tca->GetEntriesFast();++it){
1022       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1023       if(!part->IsPhysicalPrimary())continue;
1024       if(type == kTrackAODMCAll){
1025         list->Add(part);
1026         iCount++;
1027       }
1028       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1029         if(part->Charge()==0)continue;
1030         if(kTrackAODMCCharged){
1031           list->Add(part);
1032         }
1033         else {
1034           if(TMath::Abs(part->Eta())>0.9)continue;
1035           list->Add(part);
1036         }
1037         iCount++;
1038       }
1039     }
1040   }// AODMCparticle
1041   list->Sort();
1042   return iCount;
1043
1044 }