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