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