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