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