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 "AliAODJetEventBackground.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliJetKineReaderHeader.h"
  58 #include "AliGenCocktailEventHeader.h"
  59 #include "AliInputEventHandler.h"
  60
  61
  62 #include "AliAnalysisHelperJetTasks.h"
  63
  64 ClassImp(AliAnalysisTaskJetSpectrum2)
  65
  66 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
  67                                                             fJetHeaderRec(0x0),
  68                                                             fJetHeaderGen(0x0),
  69                                                             fAOD(0x0),
  70                                                             fhnCorrelation(0x0),
  71   fhnCorrelationPhiZRec(0x0),
  72   f1PtScale(0x0),
  73                                                             fBranchRec("jets"),
  74                                                             fBranchGen(""),
  75                                                             fBranchBkg(""),
  76                                                             fUseAODJetInput(kFALSE),
  77                                                             fUseAODTrackInput(kFALSE),
  78                                                             fUseAODMCInput(kFALSE),
  79                                                             fUseGlobalSelection(kFALSE),
  80                                                             fUseExternalWeightOnly(kFALSE),
  81                                                             fLimitGenJetEta(kFALSE),
  82                                                             fBkgSubtraction(kFALSE),
  83                                                             fFilterMask(0),
  84   fEventSelectionMask(0),
  85                                                             fAnalysisType(0),
  86   fTrackTypeRec(kTrackUndef),
  87   fTrackTypeGen(kTrackUndef),
  88   fAvgTrials(1),
  89   fExternalWeight(1),
  90                                                             fRecEtaWindow(0.5),
  91   fMinJetPt(0),
  92                                                             fDeltaPhiWindow(20./180.*TMath::Pi()),
  93   fh1Xsec(0x0),
  94   fh1Trials(0x0),
  95   fh1PtHard(0x0),
  96   fh1PtHardNoW(0x0),
  97   fh1PtHardTrials(0x0),
  98   fh1ZVtx(0x0),
  99   fh1NGenJets(0x0),
 100   fh1NRecJets(0x0),
 101   fh1PtTrackRec(0x0),
 102   fh1SumPtTrackRec(0x0),
 103   fh1SumPtTrackAreaRec(0x0),
 104   fh1TmpRho(0x0),
 105   fh1PtJetsRecIn(0x0),
 106   fh1PtJetsLeadingRecIn(0x0),
 107   fh1PtTracksRecIn(0x0),
 108   fh1PtTracksLeadingRecIn(0x0),
 109   fh1PtTracksGenIn(0x0),
 110   fh2NRecJetsPt(0x0),
 111   fh2NRecTracksPt(0x0),
 112   fh2JetsLeadingPhiEta(0x0),
 113   fh2JetsLeadingPhiPt(0x0),
 114   fh2TracksLeadingPhiEta(0x0),
 115   fh2TracksLeadingPhiPt(0x0),
 116   fh2TracksLeadingJetPhiPt(0x0),
 117   fh2JetPtJetPhi(0x0),
 118   fh2TrackPtTrackPhi(0x0),
 119   fh2RelPtFGen(0x0),
 120   fh2DijetDeltaPhiPt(0x0),
 121   fh2DijetAsymPt(0x0),
 122   fh2DijetAsymPtCut(0x0),
 123   fh2DijetDeltaPhiDeltaEta(0x0),
 124   fh2DijetPt2vsPt1(0x0),
 125   fh2DijetDifvsSum(0x0),
 126   fh1DijetMinv(0x0),
 127   fh1DijetMinvCut(0x0),
 128   fh1Bkg1(0x0),
 129   fh1Bkg2(0x0),
 130   fh1Bkg3(0x0),
 131   fh1Sigma1(0x0),
 132   fh1Sigma2(0x0),
 133   fh1Sigma3(0x0),
 134   fh1Area1(0x0),
 135   fh1Area2(0x0),
 136   fh1Area3(0x0),
 137   fh1Ptjet(0x0),
 138   fh1Ptjetsub1(0x0),
 139   fh1Ptjetsub2(0x0),
 140   fh1Ptjetsub3(0x0),
 141   fh1Ptjethardest(0x0),
 142   fh1Ptjetsubhardest1(0x0),
 143   fh1Ptjetsubhardest2(0x0),
 144   fh1Ptjetsubhardest3(0x0),
 145   fh1Rhovspthardest1(0x0),
 146   fh1Rhovspthardest2(0x0),
 147   fh1Rhovspthardest3(0x0),
 148   fh1Errorvspthardest1(0x0),
 149   fh1Errorvspthardest2(0x0),
 150   fh1Errorvspthardest3(0x0),
 151   fHistList(0x0)
 152 {
 153   for(int i = 0;i < kMaxStep*2;++i){
 154     fhnJetContainer[i] = 0;
 155   }
 156   for(int i = 0;i < kMaxJets;++i){
 157     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 158
 159     fh2PhiPt[i] = 0;
 160     fh2PhiEta[i] = 0;
 161     fh2RhoPtRec[i] = 0;
 162     fh2RhoPtGen[i] = 0;
 163     fh2PsiPtGen[i] = 0;
 164     fh2PsiPtRec[i] = 0;
 165     fh2FragRec[i] = 0;
 166     fh2FragLnRec[i] = 0;
 167     fh2FragGen[i] = 0;
 168     fh2FragLnGen[i] = 0;
 169   }
 170
 171 }
 172
 173 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 174   AliAnalysisTaskSE(name),
 175   fJetHeaderRec(0x0),
 176   fJetHeaderGen(0x0),
 177   fAOD(0x0),
 178   fhnCorrelation(0x0),
 179   fhnCorrelationPhiZRec(0x0),
 180   f1PtScale(0x0),
 181   fBranchRec("jets"),
 182   fBranchGen(""),
 183   fBranchBkg(""),
 184   fUseAODJetInput(kFALSE),
 185   fUseAODTrackInput(kFALSE),
 186   fUseAODMCInput(kFALSE),
 187   fUseGlobalSelection(kFALSE),
 188   fUseExternalWeightOnly(kFALSE),
 189   fLimitGenJetEta(kFALSE),
 190   fBkgSubtraction(kFALSE),
 191   fFilterMask(0),
 192   fEventSelectionMask(0),
 193   fAnalysisType(0),
 194   fTrackTypeRec(kTrackUndef),
 195   fTrackTypeGen(kTrackUndef),
 196   fAvgTrials(1),
 197   fExternalWeight(1),
 198   fRecEtaWindow(0.5),
 199   fMinJetPt(0),
 200   fDeltaPhiWindow(20./180.*TMath::Pi()),
 201   fh1Xsec(0x0),
 202   fh1Trials(0x0),
 203   fh1PtHard(0x0),
 204   fh1PtHardNoW(0x0),
 205   fh1PtHardTrials(0x0),
 206   fh1ZVtx(0x0),
 207   fh1NGenJets(0x0),
 208   fh1NRecJets(0x0),
 209   fh1PtTrackRec(0x0),
 210   fh1SumPtTrackRec(0x0),
 211   fh1SumPtTrackAreaRec(0x0),
 212   fh1TmpRho(0x0),
 213   fh1PtJetsRecIn(0x0),
 214   fh1PtJetsLeadingRecIn(0x0),
 215   fh1PtTracksRecIn(0x0),
 216   fh1PtTracksLeadingRecIn(0x0),
 217   fh1PtTracksGenIn(0x0),
 218   fh2NRecJetsPt(0x0),
 219   fh2NRecTracksPt(0x0),
 220   fh2JetsLeadingPhiEta(0x0),
 221   fh2JetsLeadingPhiPt(0x0),
 222   fh2TracksLeadingPhiEta(0x0),
 223   fh2TracksLeadingPhiPt(0x0),
 224   fh2TracksLeadingJetPhiPt(0x0),
 225   fh2JetPtJetPhi(0x0),
 226   fh2TrackPtTrackPhi(0x0),
 227   fh2RelPtFGen(0x0),
 228   fh2DijetDeltaPhiPt(0x0),
 229   fh2DijetAsymPt(0x0),
 230   fh2DijetAsymPtCut(0x0),
 231   fh2DijetDeltaPhiDeltaEta(0x0),
 232   fh2DijetPt2vsPt1(0x0),
 233   fh2DijetDifvsSum(0x0),
 234   fh1DijetMinv(0x0),
 235   fh1DijetMinvCut(0x0),
 236   fh1Bkg1(0x0),
 237   fh1Bkg2(0x0),
 238   fh1Bkg3(0x0),
 239   fh1Sigma1(0x0),
 240   fh1Sigma2(0x0),
 241   fh1Sigma3(0x0),
 242   fh1Area1(0x0),
 243   fh1Area2(0x0),
 244   fh1Area3(0x0),
 245   fh1Ptjet(0x0),
 246   fh1Ptjetsub1(0x0),
 247   fh1Ptjetsub2(0x0),
 248   fh1Ptjetsub3(0x0),
 249   fh1Ptjethardest(0x0),
 250   fh1Ptjetsubhardest1(0x0),
 251   fh1Ptjetsubhardest2(0x0),
 252   fh1Ptjetsubhardest3(0x0),
 253   fh1Rhovspthardest1(0x0),
 254   fh1Rhovspthardest2(0x0),
 255   fh1Rhovspthardest3(0x0),
 256   fh1Errorvspthardest1(0x0),
 257   fh1Errorvspthardest2(0x0),
 258   fh1Errorvspthardest3(0x0),
 259   fHistList(0x0)
 260 {
 261
 262   for(int i = 0;i < kMaxStep*2;++i){
 263     fhnJetContainer[i] = 0;
 264   }
 265   for(int i = 0;i < kMaxJets;++i){
 266     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 267
 268     fh2PhiPt[i] = 0;
 269     fh2PhiEta[i] = 0;
 270     fh2RhoPtRec[i] = 0;
 271     fh2RhoPtGen[i] = 0;
 272     fh2PsiPtGen[i] = 0;
 273     fh2PsiPtRec[i] = 0;
 274     fh2FragRec[i] = 0;
 275     fh2FragLnRec[i] = 0;
 276     fh2FragGen[i] = 0;
 277     fh2FragLnGen[i] = 0;
 278   }
 279   DefineOutput(1, TList::Class());
 280 }
 281
 282
 283
 284 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 285 {
 286   //
 287   // Implemented Notify() to read the cross sections
 288   // and number of trials from pyxsec.root
 289   //
 290
 291   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 292   Float_t xsection = 0;
 293   Float_t ftrials  = 1;
 294
 295   fAvgTrials = 1;
 296   if(tree){
 297     TFile *curfile = tree->GetCurrentFile();
 298     if (!curfile) {
 299       Error("Notify","No current file");
 300       return kFALSE;
 301     }
 302     if(!fh1Xsec||!fh1Trials){
 303       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 304       return kFALSE;
 305     }
 306     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 307     fh1Xsec->Fill("<#sigma>",xsection);
 308     // construct a poor man average trials
 309     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 310     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 311   }
 312   return kTRUE;
 313 }
 314
 315 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 316 {
 317
 318   //
 319   // Create the output container
 320   //
 321
 322
 323   // Connect the AOD
 324
 325
 326   MakeJetContainer();
 327
 328
 329   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 330
 331   OpenFile(1);
 332   if(!fHistList)fHistList = new TList();
 333   fHistList->SetOwner(kTRUE);
 334   Bool_t oldStatus = TH1::AddDirectoryStatus();
 335   TH1::AddDirectory(kFALSE);
 336
 337   //
 338   //  Histogram
 339
 340   const Int_t nBinPt = 320;
 341   Double_t binLimitsPt[nBinPt+1];
 342   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 343     if(iPt == 0){
 344       binLimitsPt[iPt] = 0.0;
 345     }
 346     else {// 1.0
 347       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 348     }
 349   }
 350
 351   const Int_t nBinPhi = 90;
 352   Double_t binLimitsPhi[nBinPhi+1];
 353   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 354     if(iPhi==0){
 355       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 356     }
 357     else{
 358       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 359     }
 360   }
 361
 362
 363   const Int_t nBinPhi2 = 360;
 364   Double_t binLimitsPhi2[nBinPhi2+1];
 365   for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
 366     if(iPhi2==0){
 367       binLimitsPhi2[iPhi2] = 0.;
 368     }
 369     else{
 370       binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
 371     }
 372   }
 373
 374
 375
 376   const Int_t nBinEta = 40;
 377   Double_t binLimitsEta[nBinEta+1];
 378   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 379     if(iEta==0){
 380       binLimitsEta[iEta] = -2.0;
 381     }
 382     else{
 383       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 384     }
 385   }
 386
 387   const Int_t nBinFrag = 25;
 388
 389
 390   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 391   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 392
 393   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 394   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 395
 396   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 397   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 398   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 399
 400   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 401
 402   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 403   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 404
 405   fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 406   fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
 407   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);
 408
 409   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 410   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 411   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 412   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 413   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 414
 415   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 416   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 417   //
 418
 419   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 420                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 421   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 422                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 423   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 424                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 425   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 426                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 427
 428   fh2JetPtJetPhi = new TH2F("fh2JetPtJetPhi","Reconstructed jet phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
 429   fh2TrackPtTrackPhi = new TH2F("fh2TrackPtTrackPhi","Reconstructed track phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
 430
 431   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-1.2,1.2);
 432
 433   for(int ij = 0;ij < kMaxJets;++ij){
 434     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 435     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 436
 437     fh2PhiPt[ij] =  new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
 438                              nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 439
 440     fh2PhiEta[ij] =  new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
 441                               nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 442
 443     fh2RhoPtRec[ij] =  new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
 444                                 40,0.,1.,nBinPt,binLimitsPt);
 445     fh2PsiPtRec[ij] =  new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
 446                                 40,0.,2.,nBinPt,binLimitsPt);
 447
 448     fh2RhoPtGen[ij] =  new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
 449                                40,0.,2.,nBinPt,binLimitsPt);
 450     fh2PsiPtGen[ij] =  new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
 451                                 40,0.,2.,nBinPt,binLimitsPt);
 452     if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 453
 454
 455     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",
 456                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 457     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",
 458                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 459
 460     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",
 461                               nBinFrag,0.,1.,nBinPt,binLimitsPt);
 462     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",
 463                                 nBinFrag,0.,10.,nBinPt,binLimitsPt);
 464   }
 465
 466   // Dijet histograms
 467
 468   fh2DijetDeltaPhiPt       = new TH2F("fh2DijetDeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
 469   fh2DijetAsymPt            = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
 470   fh2DijetAsymPtCut         = new TH2F("fh2DijetAsymCut","Pt asymmetry after delta phi cut;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
 471   fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
 472   fh2DijetPt2vsPt1          = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
 473   fh2DijetDifvsSum          = new TH2F("fh2DijetDifvsSum","Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
 474   fh1DijetMinv               = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 475   fh1DijetMinvCut           = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 476   //background histograms
 477   if(fBkgSubtraction){
 478     fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
 479     fh1Bkg2 = new TH1F("fh1Bkg2","Background estimate 2",100,0.,10.);
 480     fh1Bkg3 = new TH1F("fh1Bkg3","Background estimate 3",100,0.,10.);
 481     fh1Sigma1 = new TH1F("fh1Sigma1","Background fluctuations 1",100,0.,10.);
 482     fh1Sigma2 = new TH1F("fh1Sigma2","Background fluctuations 2",100,0.,10.);
 483     fh1Sigma3 = new TH1F("fh1Sigma3","Background fluctuations 3",100,0.,10.);
 484     fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
 485     fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
 486     fh1Area3 = new TH1F("fh1Area3","Background mean area 3",50,0.,5.);
 487     fh1Ptjet = new TH1F("fh1Ptjet","Jet spectrum",100,0.,200.);
 488     fh1Ptjetsub1 = new TH1F("fh1Ptjetsub1","Subtracted spectrum 1",50,0.,200.);
 489     fh1Ptjetsub2 = new TH1F("fh1Ptjetsub2","Subtracted spectrum 2",50,0.,200.);
 490     fh1Ptjetsub3 = new TH1F("fh1Ptjetsub3","Subtracted spectrum 3",50,0.,200.);
 491     fh1Ptjethardest = new TH1F("fh1Ptjethardest","Hardest jet spectrum",50,0.,200.);
 492     fh1Ptjetsubhardest1 = new TH1F("fh1Pthardestsub1","Subtracted hardest jet spectrum 1",100,0.,200.);
 493     fh1Ptjetsubhardest2 = new TH1F("fh1Pthardestsub2","Subtracted hardest jet spectrum 2",100,0.,200.);
 494     fh1Ptjetsubhardest3 = new TH1F("fh1Pthardestsub3","Subtracted hardest jet spectrum 3",100,0.,200.);
 495     fh1Rhovspthardest1 = new TH2F("fh1Rhovspthardest1","Background vs pTjet 1",100,0.,200.,50,0.,5.);
 496     fh1Rhovspthardest2 = new TH2F("fh1Rhovspthardest2","Background vs pTjet 2",100,0.,200.,50,0.,5.);
 497     fh1Rhovspthardest3 = new TH2F("fh1Rhovspthardest3","Background vs pTjet 3",100,0.,200.,50,0.,5.);
 498     fh1Errorvspthardest1 = new TH2F("fhErorvspthardest1","Relative error 1",100,0.,200.,50,0.,5.);
 499     fh1Errorvspthardest2 = new TH2F("fh1Errorvspthardest2","Relative error 2",100,0.,200.,50,0.,5.);
 500     fh1Errorvspthardest3 = new TH2F("fh1Errorvspthardest3","Relative error 3",100,0.,200.,50,0.,5.);
 501   }
 502
 503
 504
 505
 506   const Int_t saveLevel = 3; // large save level more histos
 507   if(saveLevel>0){
 508     fHistList->Add(fh1Xsec);
 509     fHistList->Add(fh1Trials);
 510     fHistList->Add(fh1PtHard);
 511     fHistList->Add(fh1PtHardNoW);
 512     fHistList->Add(fh1PtHardTrials);
 513     fHistList->Add(fh1ZVtx);
 514     if(fBranchGen.Length()>0){
 515       fHistList->Add(fh1NGenJets);
 516       fHistList->Add(fh1PtTracksGenIn);
 517     }
 518     fHistList->Add(fh1PtJetsRecIn);
 519     fHistList->Add(fh1PtJetsLeadingRecIn);
 520     fHistList->Add(fh1PtTracksRecIn);
 521     fHistList->Add(fh1PtTracksLeadingRecIn);
 522     fHistList->Add(fh1NRecJets);
 523     fHistList->Add(fh1PtTrackRec);
 524     fHistList->Add(fh1SumPtTrackRec);
 525     fHistList->Add(fh1SumPtTrackAreaRec);
 526     fHistList->Add(fh2NRecJetsPt);
 527     fHistList->Add(fh2NRecTracksPt);
 528     fHistList->Add(fh2JetsLeadingPhiEta );
 529     fHistList->Add(fh2JetsLeadingPhiPt );
 530     fHistList->Add(fh2TracksLeadingPhiEta);
 531     fHistList->Add(fh2TracksLeadingPhiPt);
 532     for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 533     for(int ij = 0;ij<kMaxJets;++ij){
 534       fHistList->Add( fh1PtRecIn[ij]);
 535
 536       if(fBranchGen.Length()>0){
 537         fHistList->Add(fh1PtGenIn[ij]);
 538         fHistList->Add(fh2FragGen[ij]);
 539         fHistList->Add(fh2FragLnGen[ij]);
 540         fHistList->Add(fh2RhoPtGen[ij]);
 541         fHistList->Add(fh2PsiPtGen[ij]);
 542       }
 543       fHistList->Add( fh2PhiPt[ij]);
 544       fHistList->Add( fh2PhiEta[ij]);
 545       fHistList->Add( fh2RhoPtRec[ij]);
 546       fHistList->Add( fh2PsiPtRec[ij]);
 547       fHistList->Add( fh2FragRec[ij]);
 548       fHistList->Add( fh2FragLnRec[ij]);
 549     }
 550     fHistList->Add(fhnCorrelation);
 551     fHistList->Add(fhnCorrelationPhiZRec);
 552     fHistList->Add(fh2JetPtJetPhi);
 553     fHistList->Add(fh2TrackPtTrackPhi);
 554     fHistList->Add(fh2RelPtFGen);
 555
 556     fHistList->Add(fh2DijetDeltaPhiPt);
 557     fHistList->Add(fh2DijetAsymPt);
 558     fHistList->Add(fh2DijetAsymPtCut);
 559     fHistList->Add(fh2DijetDeltaPhiDeltaEta);
 560     fHistList->Add(fh2DijetPt2vsPt1);
 561     fHistList->Add(fh2DijetDifvsSum);
 562     fHistList->Add(fh1DijetMinv);
 563     fHistList->Add(fh1DijetMinvCut);
 564     if(fBkgSubtraction){
 565       fHistList->Add(fh1Bkg1);
 566       fHistList->Add(fh1Bkg2);
 567       fHistList->Add(fh1Bkg3);
 568       fHistList->Add(fh1Sigma1);
 569       fHistList->Add(fh1Sigma2);
 570       fHistList->Add(fh1Sigma3);
 571       fHistList->Add(fh1Area1);
 572       fHistList->Add(fh1Area2);
 573       fHistList->Add(fh1Area3);
 574       fHistList->Add(fh1Ptjet);
 575       fHistList->Add(fh1Ptjethardest);
 576       fHistList->Add(fh1Ptjetsub1);
 577       fHistList->Add(fh1Ptjetsub2);
 578       fHistList->Add(fh1Ptjetsub3);
 579       fHistList->Add(fh1Ptjetsubhardest1);
 580       fHistList->Add(fh1Ptjetsubhardest2);
 581       fHistList->Add(fh1Ptjetsubhardest3);
 582       fHistList->Add(fh1Rhovspthardest1);
 583       fHistList->Add(fh1Rhovspthardest2);
 584       fHistList->Add(fh1Rhovspthardest3);
 585       fHistList->Add(fh1Errorvspthardest1);
 586       fHistList->Add(fh1Errorvspthardest2);
 587      fHistList->Add(fh1Errorvspthardest3);
 588     }
 589   }
 590
 591   // =========== Switch on Sumw2 for all histos ===========
 592   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 593     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 594     if (h1){
 595       h1->Sumw2();
 596       continue;
 597     }
 598     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 599     if(hn)hn->Sumw2();
 600   }
 601   TH1::AddDirectory(oldStatus);
 602 }
 603
 604 void AliAnalysisTaskJetSpectrum2::Init()
 605 {
 606   //
 607   // Initialization
 608   //
 609
 610   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 611
 612 }
 613
 614 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 615 {
 616
 617   Bool_t selected = kTRUE;
 618
 619   if(fUseGlobalSelection&&fEventSelectionMask==0){
 620     selected = AliAnalysisHelperJetTasks::Selected();
 621   }
 622   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 623     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 624   }
 625
 626   if(!selected){
 627     // no selection by the service task, we continue
 628     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 629     PostData(1, fHistList);
 630     return;
 631   }
 632
 633
 634   //
 635   // Execute analysis for current event
 636   //
 637   AliESDEvent *fESD = 0;
 638   if(fUseAODJetInput){
 639     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 640     if(!fAOD){
 641       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
 642       return;
 643     }
 644     // fethc the header
 645   }
 646   else{
 647     //  assume that the AOD is in the general output...
 648     fAOD  = AODEvent();
 649     if(!fAOD){
 650       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 651       return;
 652     }
 653     fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 654   }
 655
 656
 657
 658
 659   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 660
 661
 662   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 663
 664   if(!aodH){
 665     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 666     return;
 667   }
 668
 669   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 670   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 671   if(!aodRecJets){
 672     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 673     return;
 674   }
 675
 676   if(fBkgSubtraction){
 677      AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
 678
 679
 680      if(!evBkg){
 681        Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
 682        return;
 683      }
 684
 685
 686      ///just to start: some very simple plots containing rho, sigma and area of the background.
 687
 688      Int_t nJets = aodRecJets->GetEntriesFast();
 689      Float_t pthardest=0.;
 690      if(nJets!=0){
 691
 692        Float_t bkg1=evBkg->GetBackground(0);
 693        Float_t bkg2=evBkg->GetBackground(1);
 694        Float_t bkg3=evBkg->GetBackground(2);
 695        Float_t sigma1=evBkg->GetSigma(0);
 696        Float_t sigma2=evBkg->GetSigma(1);
 697        Float_t sigma3=evBkg->GetSigma(2);
 698        Float_t area1=evBkg->GetMeanarea(0);
 699        Float_t area2=evBkg->GetMeanarea(1);
 700        Float_t area3=evBkg->GetMeanarea(2);
 701        fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
 702        fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
 703        fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
 704        fh1Sigma1->Fill(sigma1);
 705        fh1Sigma2->Fill(sigma2);
 706        fh1Sigma3->Fill(sigma3);
 707        fh1Area1->Fill(area1);
 708        fh1Area2->Fill(area2);
 709        fh1Area3->Fill(area3);
 710
 711
 712        for(Int_t k=0;k<nJets;k++){
 713          AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
 714          fh1Ptjet->Fill(jet->Pt());
 715          Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
 716          Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
 717          Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
 718          Float_t err1=sigma1*sqrt(area1);
 719          Float_t err2=sigma2*sqrt(area2);
 720          Float_t err3=sigma3*sqrt(area3);
 721          fh1Ptjetsub1->Fill(ptsub1);
 722          fh1Ptjetsub2->Fill(ptsub2);
 723          fh1Ptjetsub3->Fill(ptsub3);
 724          if(k==0) {
 725            pthardest=jet->Pt();
 726            fh1Ptjethardest->Fill(pthardest);
 727            fh1Ptjetsubhardest1->Fill(ptsub1);
 728            fh1Ptjetsubhardest2->Fill(ptsub2);
 729            fh1Ptjetsubhardest3->Fill(ptsub3);
 730            fh1Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
 731            fh1Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
 732            fh1Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
 733          }
 734        }
 735        fh1Rhovspthardest1->Fill(pthardest,bkg1);
 736        fh1Rhovspthardest2->Fill(pthardest,bkg2);
 737        fh1Rhovspthardest3->Fill(pthardest,bkg3);
 738
 739
 740
 741      }
 742   }// background subtraction
 743
 744
 745   // ==== General variables needed
 746
 747
 748   // We use statice array, not to fragment the memory
 749   AliAODJet genJets[kMaxJets];
 750   Int_t nGenJets = 0;
 751   AliAODJet recJets[kMaxJets];
 752   Int_t nRecJets = 0;
 753   ///////////////////////////
 754
 755
 756   Double_t eventW = 1;
 757   Double_t ptHard = 0;
 758   Double_t nTrials = 1; // Trials for MC trigger
 759
 760   if(fUseExternalWeightOnly){
 761     eventW = fExternalWeight;
 762   }
 763
 764   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 765   //  if(fDebug>0)aodH->SetFillAOD(kFALSE);
 766   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 767   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 768     // this is the part we only use when we have MC information
 769     AliMCEvent* mcEvent = MCEvent();
 770     //    AliStack *pStack = 0;
 771     if(!mcEvent){
 772       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 773       return;
 774     }
 775     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 776     Int_t iCount = 0;
 777     if(pythiaGenHeader){
 778       nTrials = pythiaGenHeader->Trials();
 779       ptHard  = pythiaGenHeader->GetPtHard();
 780       int iProcessType = pythiaGenHeader->ProcessType();
 781       // 11 f+f -> f+f
 782       // 12 f+barf -> f+barf
 783       // 13 f+barf -> g+g
 784       // 28 f+g -> f+g
 785       // 53 g+g -> f+barf
 786       // 68 g+g -> g+g
 787       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 788       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 789
 790       // fetch the pythia generated jets only to be used here
 791       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 792       AliAODJet pythiaGenJets[kMaxJets];
 793       for(int ip = 0;ip < nPythiaGenJets;++ip){
 794         if(iCount>=kMaxJets)continue;
 795         Float_t p[4];
 796         pythiaGenHeader->TriggerJet(ip,p);
 797         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 798
 799         if(fBranchGen.Length()==0){
 800           /*
 801           if(fLimitGenJetEta){
 802             if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 803                pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 804           }
 805           */
 806           if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
 807           // if we have MC particles and we do not read from the aod branch
 808           // use the pythia jets
 809           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 810           iCount++;
 811         }
 812       }
 813     }
 814     if(fBranchGen.Length()==0)nGenJets = iCount;
 815   }// (fAnalysisType&kMCESD)==kMCESD)
 816
 817
 818   // we simply fetch the tracks/mc particles as a list of AliVParticles
 819
 820   TList recParticles;
 821   TList genParticles;
 822
 823
 824
 825
 826   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 827   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 828   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 829   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 830
 831
 832   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 833   fh1PtHard->Fill(ptHard,eventW);
 834   fh1PtHardNoW->Fill(ptHard,1);
 835   fh1PtHardTrials->Fill(ptHard,nTrials);
 836   fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
 837
 838
 839   // If we set a second branch for the input jets fetch this
 840   if(fBranchGen.Length()>0){
 841     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 842     if(aodGenJets){
 843       Int_t iCount = 0;
 844       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 845         if(iCount>=kMaxJets)continue;
 846         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 847         if(!tmp)continue;
 848         /*
 849         if(fLimitGenJetEta){
 850           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 851              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 852         }
 853         */
 854         if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
 855         genJets[iCount] = *tmp;
 856         iCount++;
 857       }
 858       nGenJets = iCount;
 859     }
 860     else{
 861       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 862       if(fDebug>2)fAOD->Print();
 863     }
 864   }
 865
 866   fh1NGenJets->Fill(nGenJets);
 867   // We do not want to exceed the maximum jet number
 868   nGenJets = TMath::Min(nGenJets,kMaxJets);
 869
 870   // Fetch the reconstructed jets...
 871
 872   nRecJets = aodRecJets->GetEntries();
 873
 874   nRecJets = aodRecJets->GetEntries();
 875   fh1NRecJets->Fill(nRecJets);
 876
 877   // Do something with the all rec jets
 878   Int_t nRecOver = nRecJets;
 879
 880   // check that the jets are sorted
 881   Float_t ptOld = 999999;
 882   for(int ir = 0;ir < nRecJets;ir++){
 883     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
 884     Float_t tmpPt = tmp->Pt();
 885     if(tmpPt>ptOld){
 886       Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
 887     }
 888     ptOld = tmpPt;
 889   }
 890
 891
 892   if(nRecOver>0){
 893     TIterator *recIter = aodRecJets->MakeIterator();
 894     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
 895     Float_t pt = tmpRec->Pt();
 896     if(tmpRec){
 897       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 898         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 899         while(pt<ptCut&&tmpRec){
 900           nRecOver--;
 901           tmpRec = (AliAODJet*)(recIter->Next());
 902           if(tmpRec){
 903             pt = tmpRec->Pt();
 904           }
 905         }
 906         if(nRecOver<=0)break;
 907         fh2NRecJetsPt->Fill(ptCut,nRecOver);
 908       }
 909     }
 910     recIter->Reset();
 911     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
 912     Float_t phi = leading->Phi();
 913     if(phi<0)phi+=TMath::Pi()*2.;
 914     Float_t eta = leading->Eta();
 915     pt = leading->Pt();
 916     while((tmpRec = (AliAODJet*)(recIter->Next()))){
 917       Float_t tmpPt = tmpRec->Pt();
 918       fh1PtJetsRecIn->Fill(tmpPt);
 919       if(tmpRec==leading){
 920         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 921         continue;
 922       }
 923       // correlation
 924       Float_t tmpPhi =  tmpRec->Phi();
 925
 926       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 927       Float_t dPhi = phi - tmpRec->Phi();
 928       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 929       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 930       Float_t dEta = eta - tmpRec->Eta();
 931       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 932       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 933     }
 934     delete recIter;
 935   }
 936
 937   Int_t nTrackOver = recParticles.GetSize();
 938   // do the same for tracks and jets
 939   if(nTrackOver>0){
 940     TIterator *recIter = recParticles.MakeIterator();
 941     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 942     Float_t pt = tmpRec->Pt();
 943     //    Printf("Leading track p_t %3.3E",pt);
 944     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 945       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 946       while(pt<ptCut&&tmpRec){
 947         nTrackOver--;
 948         tmpRec = (AliVParticle*)(recIter->Next());
 949         if(tmpRec){
 950           pt = tmpRec->Pt();
 951         }
 952       }
 953       if(nTrackOver<=0)break;
 954       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 955     }
 956
 957     recIter->Reset();
 958     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 959     Float_t phi = leading->Phi();
 960     if(phi<0)phi+=TMath::Pi()*2.;
 961     Float_t eta = leading->Eta();
 962     pt  = leading->Pt();
 963     while((tmpRec = (AliVParticle*)(recIter->Next()))){
 964       Float_t tmpPt = tmpRec->Pt();
 965       fh1PtTracksRecIn->Fill(tmpPt);
 966       if(tmpRec==leading){
 967         fh1PtTracksLeadingRecIn->Fill(tmpPt);
 968         continue;
 969       }
 970       // correlation
 971       Float_t tmpPhi =  tmpRec->Phi();
 972
 973       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 974       Float_t dPhi = phi - tmpRec->Phi();
 975       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 976       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 977       Float_t dEta = eta - tmpRec->Eta();
 978       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 979       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 980     }
 981     delete recIter;
 982   }
 983
 984   if(genParticles.GetSize()){
 985     TIterator *genIter = genParticles.MakeIterator();
 986     AliVParticle *tmpGen = 0;
 987     while((tmpGen = (AliVParticle*)(genIter->Next()))){
 988       if(TMath::Abs(tmpGen->Eta())<0.9){
 989         Float_t tmpPt = tmpGen->Pt();
 990         fh1PtTracksGenIn->Fill(tmpPt);
 991       }
 992     }
 993     delete genIter;
 994   }
 995
 996   nRecJets = TMath::Min(nRecJets,kMaxJets);
 997
 998   Int_t iCountRec = 0;
 999   for(int ir = 0;ir < nRecJets;++ir){
1000     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1001     if(!tmp)continue;
1002     if(tmp->Pt()<fMinJetPt)continue;
1003     recJets[ir] = *tmp;
1004     iCountRec++;
1005   }
1006   nRecJets = iCountRec;
1007
1008
1009   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1010   // Relate the jets
1011   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
1012   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
1013
1014
1015   for(int i = 0;i<kMaxJets;++i){
1016     iGenIndex[i] = iRecIndex[i] = -1;
1017   }
1018
1019   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1020                                             iGenIndex,iRecIndex,fDebug);
1021   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1022
1023   if(fDebug){
1024     for(int i = 0;i<kMaxJets;++i){
1025       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1026       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1027     }
1028   }
1029
1030
1031
1032
1033   Double_t container[6];
1034   Double_t containerPhiZ[6];
1035
1036   // loop over generated jets
1037
1038   // radius; tmp, get from the jet header itself
1039   // at some point, how todeal woht FastJet on MC?
1040   Float_t radiusGen =0.4;
1041   Float_t radiusRec =0.4;
1042
1043   for(int ig = 0;ig < nGenJets;++ig){
1044     Double_t ptGen = genJets[ig].Pt();
1045     Double_t phiGen = genJets[ig].Phi();
1046     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1047     Double_t etaGen = genJets[ig].Eta();
1048
1049     container[3] = ptGen;
1050     container[4] = etaGen;
1051     container[5] = phiGen;
1052     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1053     Int_t ir = iRecIndex[ig];
1054
1055     if(TMath::Abs(etaGen)<fRecEtaWindow){
1056       fh1TmpRho->Reset();
1057
1058       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1059       fh1PtGenIn[ig]->Fill(ptGen,eventW);
1060       // fill the fragmentation function
1061       for(int it = 0;it<genParticles.GetEntries();++it){
1062         AliVParticle *part = (AliVParticle*)genParticles.At(it);
1063         Float_t deltaR = genJets[ig].DeltaR(part);
1064         fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1065         if(deltaR<radiusGen){
1066           Float_t z = part->Pt()/ptGen;
1067           Float_t lnz =  -1.*TMath::Log(z);
1068           fh2FragGen[ig]->Fill(z,ptGen,eventW);
1069           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1070         }
1071
1072       }
1073       Float_t rhoSum = 0;
1074       for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
1075         Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1076         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1077         rhoSum += rho;
1078         fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1079         fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
1080       }
1081     }
1082     if(ir>=0&&ir<nRecJets){
1083       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1084       Double_t etaRec = recJets[ir].Eta();
1085       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1086       if(TMath::Abs(etaRec)<fRecEtaWindow&&TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
1087     }
1088   }// loop over generated jets
1089
1090   Float_t sumPt = 0;
1091   for(int it = 0;it<recParticles.GetEntries();++it){
1092     AliVParticle *part = (AliVParticle*)recParticles.At(it);
1093     // fill sum pt and P_t of all paritles
1094     if(TMath::Abs(part->Eta())<0.9){
1095       Float_t pt = part->Pt();
1096       fh1PtTrackRec->Fill(pt,eventW);
1097       fh2TrackPtTrackPhi->Fill(pt,part->Phi());
1098       sumPt += pt;
1099     }
1100   }
1101   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1102   fh1SumPtTrackRec->Fill(sumPt,eventW);
1103
1104
1105   // loop over reconstructed jets
1106   for(int ir = 0;ir < nRecJets;++ir){
1107     Double_t etaRec = recJets[ir].Eta();
1108     Double_t ptRec = recJets[ir].Pt();
1109     Double_t phiRec = recJets[ir].Phi();
1110     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1111
1112
1113   // do something with dijets...
1114   if(ir==1){
1115     Double_t etaRec1 = recJets[0].Eta();
1116     Double_t ptRec1 = recJets[0].Pt();
1117     Double_t phiRec1 = recJets[0].Phi();
1118     if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
1119
1120
1121     if(TMath::Abs(etaRec1)<fRecEtaWindow
1122       &&TMath::Abs(etaRec)<fRecEtaWindow){
1123
1124       Float_t deltaPhi = phiRec1 - phiRec;
1125
1126       if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1127       if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1128       deltaPhi = TMath::Abs(deltaPhi);
1129       fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
1130       Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
1131       fh2DijetAsymPt->Fill(asym,ptRec1);
1132       fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
1133       fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
1134       fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
1135       Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
1136                          recJets[0].Px()*recJets[1].Px()-
1137                          recJets[0].Py()*recJets[1].Py()-
1138                          recJets[0].Pz()*recJets[1].Py());
1139       minv = TMath::Sqrt(minv);
1140       // with mass == 0;
1141
1142       fh1DijetMinv->Fill(minv);
1143       if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
1144         fh1DijetMinvCut->Fill(minv);
1145         fh2DijetAsymPtCut->Fill(asym,ptRec1);
1146       }
1147     }
1148   }
1149
1150
1151     container[0] = ptRec;
1152     container[1] = etaRec;
1153     container[2] = phiRec;
1154     containerPhiZ[0] = ptRec;
1155     containerPhiZ[1] = phiRec;
1156     if(ptRec>30.&&fDebug>0){
1157       // need to cast to int, otherwise the printf overwrites
1158       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1159       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1160       if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
1161       //  aodH->SetFillAOD(kTRUE);
1162       fAOD->GetHeader()->Print();
1163       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
1164       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1165         AliAODTrack *tr = fAOD->GetTrack(it);
1166         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1167         tr->Print();
1168         //      tr->Dump();
1169         if(fESD){
1170           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1171           esdTr->Print("");
1172           // esdTr->Dump();
1173         }
1174       }
1175     }
1176
1177
1178     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1179     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1180
1181     Float_t zLeading = -1;
1182     if(TMath::Abs(etaRec)<fRecEtaWindow){
1183       fh2JetPtJetPhi->Fill(ptRec,phiRec);
1184       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1185       fh1PtRecIn[ir]->Fill(ptRec,eventW);
1186       // fill the fragmentation function
1187
1188       fh1TmpRho->Reset();
1189
1190       for(int it = 0;it<recParticles.GetEntries();++it){
1191         AliVParticle *part = (AliVParticle*)recParticles.At(it);
1192         Float_t eta = part->Eta();
1193         if(TMath::Abs(eta)<0.9){
1194           Float_t phi = part->Phi();
1195           if(phi<0)phi+=TMath::Pi()*2.;
1196           Float_t dPhi = phi - phiRec;
1197           Float_t dEta = eta - etaRec;
1198           if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1199           if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1200           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1201           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1202         }
1203
1204         Float_t deltaR = recJets[ir].DeltaR(part);
1205         fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1206
1207
1208         if(deltaR<radiusRec){
1209           Float_t z = part->Pt()/ptRec;
1210           if(z>zLeading)zLeading=z;
1211           Float_t lnz =  -1.*TMath::Log(z);
1212           fh2FragRec[ir]->Fill(z,ptRec,eventW);
1213           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1214         }
1215       }
1216       // fill the jet shapes
1217       Float_t rhoSum = 0;
1218       for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1219         Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1220         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1221         rhoSum += rho;
1222         fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1223         fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1224       }
1225     }
1226
1227
1228     containerPhiZ[2] = zLeading;
1229
1230     // Fill Correlation
1231     Int_t ig = iGenIndex[ir];
1232     if(ig>=0 && ig<nGenJets){
1233       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1234       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1235       Double_t ptGen  = genJets[ig].Pt();
1236       Double_t phiGen = genJets[ig].Phi();
1237       if(phiGen<0)phiGen+=TMath::Pi()*2.;
1238       Double_t etaGen = genJets[ig].Eta();
1239
1240       container[3] = ptGen;
1241       container[4] = etaGen;
1242       container[5] = phiGen;
1243       containerPhiZ[3] = ptGen;
1244       //
1245       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1246       //
1247
1248       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1249       if(TMath::Abs(etaRec)<fRecEtaWindow){
1250         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1251         fhnCorrelation->Fill(container);
1252         if(ptGen>0){
1253           Float_t delta = (ptRec-ptGen)/ptGen;
1254           fh2RelPtFGen->Fill(ptGen,delta,eventW);
1255         }
1256         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1257
1258       }// if etarec in window
1259
1260     }
1261     else{
1262       containerPhiZ[3] = 0;
1263       if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1264     }
1265   }// loop over reconstructed jets
1266
1267
1268   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1269   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1270   PostData(1, fHistList);
1271 }
1272
1273
1274 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1275   //
1276   // Create the particle container for the correction framework manager and
1277   // link it
1278   //
1279   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1280   const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1281   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1282   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1283   const Double_t kZmin = 0., kZmax = 1;
1284
1285   // can we neglect migration in eta and phi?
1286   // phi should be no problem since we cover full phi and are phi symmetric
1287   // eta migration is more difficult  due to needed acceptance correction
1288   // in limited eta range
1289
1290   //arrays for the number of bins in each dimension
1291   Int_t iBin[kNvar];
1292   iBin[0] = 320; //bins in pt
1293   iBin[1] =  1; //bins in eta
1294   iBin[2] = 1; // bins in phi
1295
1296   //arrays for lower bounds :
1297   Double_t* binEdges[kNvar];
1298   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1299     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1300
1301   //values for bin lower bounds
1302   //  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);
1303   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1304   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1305   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1306
1307
1308   for(int i = 0;i < kMaxStep*2;++i){
1309     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1310     for (int k=0; k<kNvar; k++) {
1311       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1312     }
1313   }
1314   //create correlation matrix for unfolding
1315   Int_t thnDim[2*kNvar];
1316   for (int k=0; k<kNvar; k++) {
1317     //first half  : reconstructed
1318     //second half : MC
1319     thnDim[k]      = iBin[k];
1320     thnDim[k+kNvar] = iBin[k];
1321   }
1322
1323   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1324   for (int k=0; k<kNvar; k++) {
1325     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1326     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1327   }
1328   fhnCorrelation->Sumw2();
1329
1330   // for second correlation histogram
1331
1332
1333   const Int_t kNvarPhiZ   = 4;
1334   //arrays for the number of bins in each dimension
1335   Int_t iBinPhiZ[kNvarPhiZ];
1336   iBinPhiZ[0] = 80; //bins in pt
1337   iBinPhiZ[1] = 72; //bins in phi
1338   iBinPhiZ[2] = 20; // bins in Z
1339   iBinPhiZ[3] = 80; //bins in ptgen
1340
1341   //arrays for lower bounds :
1342   Double_t* binEdgesPhiZ[kNvarPhiZ];
1343   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1344     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1345
1346   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1347   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1348   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1349   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1350
1351   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1352   for (int k=0; k<kNvarPhiZ; k++) {
1353     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1354   }
1355   fhnCorrelationPhiZRec->Sumw2();
1356
1357
1358   // Add a histogram for Fake jets
1359
1360   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1361     delete [] binEdges[ivar];
1362
1363   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1364     delete [] binEdgesPhiZ[ivar];
1365
1366 }
1367
1368 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1369 {
1370 // Terminate analysis
1371 //
1372     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1373 }
1374
1375
1376 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1377
1378   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1379
1380   Int_t iCount = 0;
1381   if(type==kTrackAOD){
1382     AliAODEvent *aod = 0;
1383     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1384     else aod = AODEvent();
1385     if(!aod){
1386       return iCount;
1387     }
1388     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1389       AliAODTrack *tr = aod->GetTrack(it);
1390       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1391       if(TMath::Abs(tr->Eta())>0.9)continue;
1392       if(fDebug>0){
1393         if(tr->Pt()>20){
1394           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1395           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1396           tr->Print();
1397           //    tr->Dump();
1398           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1399           if(fESD){
1400             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1401             esdTr->Print("");
1402             // esdTr->Dump();
1403           }
1404         }
1405       }
1406       list->Add(tr);
1407       iCount++;
1408     }
1409   }
1410   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1411     AliMCEvent* mcEvent = MCEvent();
1412     if(!mcEvent)return iCount;
1413     // we want to have alivpartilces so use get track
1414     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1415       if(!mcEvent->IsPhysicalPrimary(it))continue;
1416       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1417       if(type == kTrackKineAll){
1418         list->Add(part);
1419         iCount++;
1420       }
1421       else if(type == kTrackKineCharged){
1422         if(part->Particle()->GetPDG()->Charge()==0)continue;
1423         list->Add(part);
1424         iCount++;
1425       }
1426     }
1427   }
1428   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1429     AliAODEvent *aod = 0;
1430     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1431     else aod = AODEvent();
1432     if(!aod)return iCount;
1433     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1434     if(!tca)return iCount;
1435     for(int it = 0;it < tca->GetEntriesFast();++it){
1436       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1437       if(!part->IsPhysicalPrimary())continue;
1438       if(type == kTrackAODMCAll){
1439         list->Add(part);
1440         iCount++;
1441       }
1442       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1443         if(part->Charge()==0)continue;
1444         if(kTrackAODMCCharged){
1445           list->Add(part);
1446         }
1447         else {
1448           if(TMath::Abs(part->Eta())>0.9)continue;
1449           list->Add(part);
1450         }
1451         iCount++;
1452       }
1453     }
1454   }// AODMCparticle
1455   list->Sort();
1456   return iCount;
1457
1458 }