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       }
 456       fHistList->Add( fh2PhiPt[ij]);
 457       fHistList->Add( fh2PhiEta[ij]);
 458       fHistList->Add( fh2RhoPtRec[ij]);
 459       fHistList->Add( fh2PsiPtRec[ij]);
 460       fHistList->Add( fh2FragRec[ij]);
 461       fHistList->Add( fh2FragLnRec[ij]);
 462     }
 463     fHistList->Add(fhnCorrelation);
 464     fHistList->Add(fhnCorrelationPhiZRec);
 465     fHistList->Add(fh2JetPtJetPhi);
 466     fHistList->Add(fh2TrackPtTrackPhi);
 467
 468     fHistList->Add(fh2DijetDeltaPhiPt);
 469     fHistList->Add(fh2DijetAsymPt);
 470     fHistList->Add(fh2DijetAsymPtCut);
 471     fHistList->Add(fh2DijetDeltaPhiDeltaEta);
 472     fHistList->Add(fh2DijetPt2vsPt1);
 473     fHistList->Add(fh2DijetDifvsSum);
 474     fHistList->Add(fh1DijetMinv);
 475     fHistList->Add(fh1DijetMinvCut);
 476   }
 477
 478   // =========== Switch on Sumw2 for all histos ===========
 479   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 480     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 481     if (h1){
 482       h1->Sumw2();
 483       continue;
 484     }
 485     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 486     if(hn)hn->Sumw2();
 487   }
 488   TH1::AddDirectory(oldStatus);
 489 }
 490
 491 void AliAnalysisTaskJetSpectrum2::Init()
 492 {
 493   //
 494   // Initialization
 495   //
 496
 497   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 498
 499 }
 500
 501 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 502 {
 503
 504   if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
 505     // no selection by the service task, we continue
 506     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 507     PostData(1, fHistList);
 508     return;
 509   }
 510
 511   //
 512   // Execute analysis for current event
 513   //
 514   AliESDEvent *fESD = 0;
 515   if(fUseAODJetInput){
 516     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 517     if(!fAOD){
 518       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
 519       return;
 520     }
 521     // fethc the header
 522   }
 523   else{
 524     //  assume that the AOD is in the general output...
 525     fAOD  = AODEvent();
 526     if(!fAOD){
 527       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 528       return;
 529     }
 530     if(fDebug>0){
 531       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 532     }
 533   }
 534
 535
 536
 537
 538   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 539
 540
 541   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 542
 543   if(!aodH){
 544     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 545     return;
 546   }
 547
 548   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 549   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 550   if(!aodRecJets){
 551     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 552     return;
 553   }
 554
 555   // ==== General variables needed
 556
 557
 558   // We use statice array, not to fragment the memory
 559   AliAODJet genJets[kMaxJets];
 560   Int_t nGenJets = 0;
 561   AliAODJet recJets[kMaxJets];
 562   Int_t nRecJets = 0;
 563   ///////////////////////////
 564
 565
 566   Double_t eventW = 1;
 567   Double_t ptHard = 0;
 568   Double_t nTrials = 1; // Trials for MC trigger
 569
 570   if(fUseExternalWeightOnly){
 571     eventW = fExternalWeight;
 572   }
 573
 574   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 575   //  if(fDebug>0)aodH->SetFillAOD(kFALSE);
 576   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 577   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 578     // this is the part we only use when we have MC information
 579     AliMCEvent* mcEvent = MCEvent();
 580     //    AliStack *pStack = 0;
 581     if(!mcEvent){
 582       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 583       return;
 584     }
 585     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 586     Int_t iCount = 0;
 587     if(pythiaGenHeader){
 588       nTrials = pythiaGenHeader->Trials();
 589       ptHard  = pythiaGenHeader->GetPtHard();
 590       int iProcessType = pythiaGenHeader->ProcessType();
 591       // 11 f+f -> f+f
 592       // 12 f+barf -> f+barf
 593       // 13 f+barf -> g+g
 594       // 28 f+g -> f+g
 595       // 53 g+g -> f+barf
 596       // 68 g+g -> g+g
 597       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 598       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 599
 600       // fetch the pythia generated jets only to be used here
 601       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 602       AliAODJet pythiaGenJets[kMaxJets];
 603       for(int ip = 0;ip < nPythiaGenJets;++ip){
 604         if(iCount>=kMaxJets)continue;
 605         Float_t p[4];
 606         pythiaGenHeader->TriggerJet(ip,p);
 607         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 608
 609         if(fBranchGen.Length()==0){
 610           /*
 611           if(fLimitGenJetEta){
 612             if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 613                pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 614           }
 615           */
 616           if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
 617           // if we have MC particles and we do not read from the aod branch
 618           // use the pythia jets
 619           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 620           iCount++;
 621         }
 622       }
 623     }
 624     if(fBranchGen.Length()==0)nGenJets = iCount;
 625   }// (fAnalysisType&kMCESD)==kMCESD)
 626
 627
 628   // we simply fetch the tracks/mc particles as a list of AliVParticles
 629
 630   TList recParticles;
 631   TList genParticles;
 632
 633
 634
 635
 636   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 637   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 638   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 639   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 640
 641
 642   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 643   fh1PtHard->Fill(ptHard,eventW);
 644   fh1PtHardNoW->Fill(ptHard,1);
 645   fh1PtHardTrials->Fill(ptHard,nTrials);
 646
 647   // If we set a second branch for the input jets fetch this
 648   if(fBranchGen.Length()>0){
 649     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 650     if(aodGenJets){
 651       Int_t iCount = 0;
 652       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 653         if(iCount>=kMaxJets)continue;
 654         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 655         if(!tmp)continue;
 656         /*
 657         if(fLimitGenJetEta){
 658           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 659              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 660         }
 661         */
 662         if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
 663         genJets[iCount] = *tmp;
 664         iCount++;
 665       }
 666       nGenJets = iCount;
 667     }
 668     else{
 669       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 670       if(fDebug>2)fAOD->Print();
 671     }
 672   }
 673
 674   fh1NGenJets->Fill(nGenJets);
 675   // We do not want to exceed the maximum jet number
 676   nGenJets = TMath::Min(nGenJets,kMaxJets);
 677
 678   // Fetch the reconstructed jets...
 679
 680   nRecJets = aodRecJets->GetEntries();
 681
 682   nRecJets = aodRecJets->GetEntries();
 683   fh1NRecJets->Fill(nRecJets);
 684
 685   // Do something with the all rec jets
 686   Int_t nRecOver = nRecJets;
 687
 688   // check that the jets are sorted
 689   Float_t ptOld = 999999;
 690   for(int ir = 0;ir < nRecJets;ir++){
 691     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
 692     Float_t tmpPt = tmp->Pt();
 693     if(tmpPt>ptOld){
 694       Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
 695     }
 696     ptOld = tmpPt;
 697   }
 698
 699
 700   if(nRecOver>0){
 701     TIterator *recIter = aodRecJets->MakeIterator();
 702     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
 703     Float_t pt = tmpRec->Pt();
 704     if(tmpRec){
 705       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 706         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 707         while(pt<ptCut&&tmpRec){
 708           nRecOver--;
 709           tmpRec = (AliAODJet*)(recIter->Next());
 710           if(tmpRec){
 711             pt = tmpRec->Pt();
 712           }
 713         }
 714         if(nRecOver<=0)break;
 715         fh2NRecJetsPt->Fill(ptCut,nRecOver);
 716       }
 717     }
 718     recIter->Reset();
 719     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
 720     Float_t phi = leading->Phi();
 721     if(phi<0)phi+=TMath::Pi()*2.;
 722     Float_t eta = leading->Eta();
 723     pt = leading->Pt();
 724      while((tmpRec = (AliAODJet*)(recIter->Next()))){
 725       Float_t tmpPt = tmpRec->Pt();
 726       fh1PtJetsRecIn->Fill(tmpPt);
 727       if(tmpRec==leading){
 728         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 729         continue;
 730       }
 731       // correlation
 732       Float_t tmpPhi =  tmpRec->Phi();
 733
 734       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 735       Float_t dPhi = phi - tmpRec->Phi();
 736       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 737       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 738       Float_t dEta = eta - tmpRec->Eta();
 739       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 740       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 741     }
 742     delete recIter;
 743   }
 744
 745   Int_t nTrackOver = recParticles.GetSize();
 746   // do the same for tracks and jets
 747   if(nTrackOver>0){
 748     TIterator *recIter = recParticles.MakeIterator();
 749     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 750     Float_t pt = tmpRec->Pt();
 751     //    Printf("Leading track p_t %3.3E",pt);
 752     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 753       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 754       while(pt<ptCut&&tmpRec){
 755         nTrackOver--;
 756         tmpRec = (AliVParticle*)(recIter->Next());
 757         if(tmpRec){
 758           pt = tmpRec->Pt();
 759         }
 760       }
 761       if(nTrackOver<=0)break;
 762       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 763     }
 764
 765     recIter->Reset();
 766     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 767     Float_t phi = leading->Phi();
 768     if(phi<0)phi+=TMath::Pi()*2.;
 769     Float_t eta = leading->Eta();
 770     pt  = leading->Pt();
 771     while((tmpRec = (AliVParticle*)(recIter->Next()))){
 772       Float_t tmpPt = tmpRec->Pt();
 773       fh1PtTracksRecIn->Fill(tmpPt);
 774       if(tmpRec==leading){
 775         fh1PtTracksLeadingRecIn->Fill(tmpPt);
 776         continue;
 777       }
 778       // correlation
 779       Float_t tmpPhi =  tmpRec->Phi();
 780
 781       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 782       Float_t dPhi = phi - tmpRec->Phi();
 783       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 784       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 785       Float_t dEta = eta - tmpRec->Eta();
 786       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 787       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 788     }
 789     delete recIter;
 790   }
 791
 792   if(genParticles.GetSize()){
 793     TIterator *genIter = genParticles.MakeIterator();
 794     AliVParticle *tmpGen = 0;
 795     while((tmpGen = (AliVParticle*)(genIter->Next()))){
 796       if(TMath::Abs(tmpGen->Eta())<0.9){
 797         Float_t tmpPt = tmpGen->Pt();
 798         fh1PtTracksGenIn->Fill(tmpPt);
 799       }
 800     }
 801     delete genIter;
 802   }
 803
 804   nRecJets = TMath::Min(nRecJets,kMaxJets);
 805
 806   Int_t iCountRec = 0;
 807   for(int ir = 0;ir < nRecJets;++ir){
 808     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 809     if(!tmp)continue;
 810     if(tmp->Pt()<fMinJetPt)continue;
 811     recJets[ir] = *tmp;
 812     iCountRec++;
 813   }
 814   nRecJets = iCountRec;
 815
 816
 817   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 818   // Relate the jets
 819   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 820   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 821
 822
 823   for(int i = 0;i<kMaxJets;++i){
 824     iGenIndex[i] = iRecIndex[i] = -1;
 825   }
 826
 827   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 828                                             iGenIndex,iRecIndex,fDebug);
 829   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 830
 831   if(fDebug){
 832     for(int i = 0;i<kMaxJets;++i){
 833       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 834       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 835     }
 836   }
 837
 838
 839
 840
 841   Double_t container[6];
 842   Double_t containerPhiZ[6];
 843
 844   // loop over generated jets
 845
 846   // radius; tmp, get from the jet header itself
 847   // at some point, how todeal woht FastJet on MC?
 848   Float_t radiusGen =0.4;
 849   Float_t radiusRec =0.4;
 850
 851   for(int ig = 0;ig < nGenJets;++ig){
 852     Double_t ptGen = genJets[ig].Pt();
 853     Double_t phiGen = genJets[ig].Phi();
 854     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 855     Double_t etaGen = genJets[ig].Eta();
 856
 857     container[3] = ptGen;
 858     container[4] = etaGen;
 859     container[5] = phiGen;
 860     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
 861     Int_t ir = iRecIndex[ig];
 862
 863     if(TMath::Abs(etaGen)<fRecEtaWindow){
 864       fh1TmpRho->Reset();
 865
 866       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
 867       fh1PtGenIn[ig]->Fill(ptGen,eventW);
 868       // fill the fragmentation function
 869       for(int it = 0;it<genParticles.GetEntries();++it){
 870         AliVParticle *part = (AliVParticle*)genParticles.At(it);
 871         Float_t deltaR = genJets[ig].DeltaR(part);
 872         fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
 873         if(deltaR<radiusGen){
 874           Float_t z = part->Pt()/ptGen;
 875           Float_t lnz =  -1.*TMath::Log(z);
 876           fh2FragGen[ig]->Fill(z,ptGen,eventW);
 877           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
 878         }
 879
 880       }
 881       Float_t rhoSum = 0;
 882       for(int ibx = 1;ibx<fh2RhoPtGen[ir]->GetNbinsX();ibx++){
 883         Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
 884         Float_t rho = fh1TmpRho->GetBinContent(ibx);
 885         rhoSum += rho;
 886         fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
 887         fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
 888       }
 889     }
 890     if(ir>=0&&ir<nRecJets){
 891       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
 892       Double_t etaRec = recJets[ir].Eta();
 893       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
 894     }
 895   }// loop over generated jets
 896
 897   Float_t sumPt = 0;
 898   for(int it = 0;it<recParticles.GetEntries();++it){
 899     AliVParticle *part = (AliVParticle*)recParticles.At(it);
 900     // fill sum pt and P_t of all paritles
 901     if(TMath::Abs(part->Eta())<0.9){
 902       Float_t pt = part->Pt();
 903       fh1PtTrackRec->Fill(pt,eventW);
 904       fh2TrackPtTrackPhi->Fill(pt,part->Phi());
 905       sumPt += pt;
 906     }
 907   }
 908   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
 909   fh1SumPtTrackRec->Fill(sumPt,eventW);
 910
 911
 912   // loop over reconstructed jets
 913   for(int ir = 0;ir < nRecJets;++ir){
 914     Double_t etaRec = recJets[ir].Eta();
 915     Double_t ptRec = recJets[ir].Pt();
 916     Double_t phiRec = recJets[ir].Phi();
 917     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 918
 919
 920   // do something with dijets...
 921   if(ir==1){
 922     Double_t etaRec1 = recJets[0].Eta();
 923     Double_t ptRec1 = recJets[0].Pt();
 924     Double_t phiRec1 = recJets[0].Phi();
 925     if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
 926
 927
 928     if(TMath::Abs(etaRec1)<fRecEtaWindow
 929       &&TMath::Abs(etaRec)<fRecEtaWindow){
 930
 931       Float_t deltaPhi = phiRec1 - phiRec;
 932
 933       if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
 934       if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
 935       deltaPhi = TMath::Abs(deltaPhi);
 936       fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
 937       Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
 938       fh2DijetAsymPt->Fill(asym,ptRec1);
 939       fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
 940       fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
 941       fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
 942       Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
 943                          recJets[0].Px()*recJets[1].Px()-
 944                          recJets[0].Py()*recJets[1].Py()-
 945                          recJets[0].Pz()*recJets[1].Py());
 946       minv = TMath::Sqrt(minv);
 947       // with mass == 0;
 948
 949       fh1DijetMinv->Fill(minv);
 950       if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
 951         fh1DijetMinvCut->Fill(minv);
 952         fh2DijetAsymPtCut->Fill(asym,ptRec1);
 953       }
 954     }
 955   }
 956
 957
 958     container[0] = ptRec;
 959     container[1] = etaRec;
 960     container[2] = phiRec;
 961     containerPhiZ[0] = ptRec;
 962     containerPhiZ[1] = phiRec;
 963     if(ptRec>30.&&fDebug>0){
 964       // need to cast to int, otherwise the printf overwrites
 965       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
 966       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
 967       if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
 968       //  aodH->SetFillAOD(kTRUE);
 969       fAOD->GetHeader()->Print();
 970       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
 971       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
 972         AliAODTrack *tr = fAOD->GetTrack(it);
 973         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 974         tr->Print();
 975         //      tr->Dump();
 976         if(fESD){
 977           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
 978           esdTr->Print("");
 979           // esdTr->Dump();
 980         }
 981       }
 982     }
 983
 984
 985     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
 986     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 987
 988     Float_t zLeading = -1;
 989     if(TMath::Abs(etaRec)<fRecEtaWindow){
 990       fh2JetPtJetPhi->Fill(ptRec,phiRec);
 991       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
 992       fh1PtRecIn[ir]->Fill(ptRec,eventW);
 993       // fill the fragmentation function
 994
 995       fh1TmpRho->Reset();
 996
 997       for(int it = 0;it<recParticles.GetEntries();++it){
 998         AliVParticle *part = (AliVParticle*)recParticles.At(it);
 999         Float_t eta = part->Eta();
1000         if(TMath::Abs(eta)<0.9){
1001           Float_t phi = part->Phi();
1002           if(phi<0)phi+=TMath::Pi()*2.;
1003           Float_t dPhi = phi - phiRec;
1004           Float_t dEta = eta - etaRec;
1005           if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1006           if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1007           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1008           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1009         }
1010
1011         Float_t deltaR = recJets[ir].DeltaR(part);
1012         fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1013
1014
1015         if(deltaR<radiusRec){
1016           Float_t z = part->Pt()/ptRec;
1017           if(z>zLeading)zLeading=z;
1018           Float_t lnz =  -1.*TMath::Log(z);
1019           fh2FragRec[ir]->Fill(z,ptRec,eventW);
1020           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1021         }
1022       }
1023       // fill the jet shapes
1024       Float_t rhoSum = 0;
1025       for(int ibx = 1;ibx<fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1026         Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1027         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1028         rhoSum += rho;
1029         fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1030         fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1031       }
1032     }
1033
1034
1035     containerPhiZ[2] = zLeading;
1036
1037     // Fill Correlation
1038     Int_t ig = iGenIndex[ir];
1039     if(ig>=0 && ig<nGenJets){
1040       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1041       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1042       Double_t ptGen  = genJets[ig].Pt();
1043       Double_t phiGen = genJets[ig].Phi();
1044       if(phiGen<0)phiGen+=TMath::Pi()*2.;
1045       Double_t etaGen = genJets[ig].Eta();
1046
1047       container[3] = ptGen;
1048       container[4] = etaGen;
1049       container[5] = phiGen;
1050       containerPhiZ[3] = ptGen;
1051       //
1052       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1053       //
1054
1055       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1056       if(TMath::Abs(etaRec)<fRecEtaWindow){
1057         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1058         fhnCorrelation->Fill(container);
1059         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1060
1061       }// if etarec in window
1062
1063     }
1064     else{
1065       containerPhiZ[3] = 0;
1066       if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1067     }
1068   }// loop over reconstructed jets
1069
1070
1071   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1072   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1073   PostData(1, fHistList);
1074 }
1075
1076
1077 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1078   //
1079   // Create the particle container for the correction framework manager and
1080   // link it
1081   //
1082   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1083   const Double_t kPtmin = 0.0, kPtmax = 160.; // we do not want to have empty bins at the beginning...
1084   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1085   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1086   const Double_t kZmin = 0., kZmax = 1;
1087
1088   // can we neglect migration in eta and phi?
1089   // phi should be no problem since we cover full phi and are phi symmetric
1090   // eta migration is more difficult  due to needed acceptance correction
1091   // in limited eta range
1092
1093   //arrays for the number of bins in each dimension
1094   Int_t iBin[kNvar];
1095   iBin[0] = 160; //bins in pt
1096   iBin[1] =  1; //bins in eta
1097   iBin[2] = 1; // bins in phi
1098
1099   //arrays for lower bounds :
1100   Double_t* binEdges[kNvar];
1101   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1102     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1103
1104   //values for bin lower bounds
1105   //  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);
1106   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1107   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1108   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1109
1110
1111   for(int i = 0;i < kMaxStep*2;++i){
1112     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
1113     for (int k=0; k<kNvar; k++) {
1114       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1115     }
1116   }
1117   //create correlation matrix for unfolding
1118   Int_t thnDim[2*kNvar];
1119   for (int k=0; k<kNvar; k++) {
1120     //first half  : reconstructed
1121     //second half : MC
1122     thnDim[k]      = iBin[k];
1123     thnDim[k+kNvar] = iBin[k];
1124   }
1125
1126   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1127   for (int k=0; k<kNvar; k++) {
1128     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1129     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1130   }
1131   fhnCorrelation->Sumw2();
1132
1133   // for second correlation histogram
1134
1135
1136   const Int_t kNvarPhiZ   = 4;
1137   //arrays for the number of bins in each dimension
1138   Int_t iBinPhiZ[kNvarPhiZ];
1139   iBinPhiZ[0] = 80; //bins in pt
1140   iBinPhiZ[1] = 72; //bins in phi
1141   iBinPhiZ[2] = 20; // bins in Z
1142   iBinPhiZ[3] = 80; //bins in ptgen
1143
1144   //arrays for lower bounds :
1145   Double_t* binEdgesPhiZ[kNvarPhiZ];
1146   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1147     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1148
1149   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1150   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1151   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1152   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1153
1154   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1155   for (int k=0; k<kNvarPhiZ; k++) {
1156     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1157   }
1158   fhnCorrelationPhiZRec->Sumw2();
1159
1160
1161   // Add a histogram for Fake jets
1162
1163   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1164     delete [] binEdges[ivar];
1165
1166   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1167     delete [] binEdgesPhiZ[ivar];
1168
1169 }
1170
1171 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1172 {
1173 // Terminate analysis
1174 //
1175     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1176 }
1177
1178
1179 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1180
1181   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1182
1183   Int_t iCount = 0;
1184   if(type==kTrackAOD){
1185     AliAODEvent *aod = 0;
1186     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1187     else aod = AODEvent();
1188     if(!aod){
1189       return iCount;
1190     }
1191     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1192       AliAODTrack *tr = aod->GetTrack(it);
1193       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1194       if(TMath::Abs(tr->Eta())>0.9)continue;
1195       if(fDebug>0){
1196         if(tr->Pt()>20){
1197           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1198           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
1199           tr->Print();
1200           //    tr->Dump();
1201           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1202           if(fESD){
1203             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1204             esdTr->Print("");
1205             // esdTr->Dump();
1206           }
1207         }
1208       }
1209       list->Add(tr);
1210       iCount++;
1211     }
1212   }
1213   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1214     AliMCEvent* mcEvent = MCEvent();
1215     if(!mcEvent)return iCount;
1216     // we want to have alivpartilces so use get track
1217     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1218       if(!mcEvent->IsPhysicalPrimary(it))continue;
1219       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1220       if(type == kTrackKineAll){
1221         list->Add(part);
1222         iCount++;
1223       }
1224       else if(type == kTrackKineCharged){
1225         if(part->Particle()->GetPDG()->Charge()==0)continue;
1226         list->Add(part);
1227         iCount++;
1228       }
1229     }
1230   }
1231   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1232     AliAODEvent *aod = 0;
1233     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1234     else aod = AODEvent();
1235     if(!aod)return iCount;
1236     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1237     if(!tca)return iCount;
1238     for(int it = 0;it < tca->GetEntriesFast();++it){
1239       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1240       if(!part->IsPhysicalPrimary())continue;
1241       if(type == kTrackAODMCAll){
1242         list->Add(part);
1243         iCount++;
1244       }
1245       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1246         if(part->Charge()==0)continue;
1247         if(kTrackAODMCCharged){
1248           list->Add(part);
1249         }
1250         else {
1251           if(TMath::Abs(part->Eta())>0.9)continue;
1252           list->Add(part);
1253         }
1254         iCount++;
1255       }
1256     }
1257   }// AODMCparticle
1258   list->Sort();
1259   return iCount;
1260
1261 }