PWG4/JetTasks/AliAnalysisTaskJetSpectrum.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 <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TH3F.h>
  32 #include <TProfile.h>
  33 #include <TList.h>
  34 #include <TLorentzVector.h>
  35 #include <TClonesArray.h>
  36 #include  "TDatabasePDG.h"
  37
  38 #include "AliAnalysisTaskJetSpectrum.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliJetFinder.h"
  41 #include "AliJetHeader.h"
  42 #include "AliJetReader.h"
  43 #include "AliJetReaderHeader.h"
  44 #include "AliUA1JetHeaderV1.h"
  45 #include "AliJet.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODTrack.h"
  50 #include "AliAODJet.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliGenPythiaEventHeader.h"
  55 #include "AliJetKineReaderHeader.h"
  56 #include "AliGenCocktailEventHeader.h"
  57 #include "AliInputEventHandler.h"
  58
  59
  60 #include "AliAnalysisHelperJetTasks.h"
  61
  62 ClassImp(AliAnalysisTaskJetSpectrum)
  63
  64 const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
  65
  66 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
  67   fJetHeaderRec(0x0),
  68   fJetHeaderGen(0x0),
  69   fAOD(0x0),
  70   fBranchRec("jets"),
  71   fBranchGen(""),
  72   fUseAODInput(kFALSE),
  73   fUseExternalWeightOnly(kFALSE),
  74   fLimitGenJetEta(kFALSE),
  75   fAnalysisType(0),
  76   fExternalWeight(1),
  77   fRecEtaWindow(0.5),
  78   fh1Xsec(0x0),
  79   fh1Trials(0x0),
  80   fh1PtHard(0x0),
  81   fh1PtHardNoW(0x0),
  82   fh1PtHardTrials(0x0),
  83   fh1NGenJets(0x0),
  84   fh1NRecJets(0x0),
  85   fHistList(0x0) ,
  86   ////////////////
  87   fh1JetMultiplicity(0x0) ,
  88   fh2ERecZRec(0x0) ,
  89   fh2EGenZGen(0x0) ,
  90   fh2Efficiency(0x0) ,
  91   fh3EGenERecN(0x0)
  92   ////////////////
  93 {
  94   // Default constructor
  95     for(int ij  = 0;ij<kMaxJets;++ij){
  96       fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
  97       fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =  fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
  98       fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
  99     }
 100     for(int ic = 0;ic < kMaxCorrelation;ic++){
 101       fhnCorrelation[ic] = 0;
 102     }
 103
 104 }
 105
 106 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
 107   AliAnalysisTaskSE(name),
 108   fJetHeaderRec(0x0),
 109   fJetHeaderGen(0x0),
 110   fAOD(0x0),
 111   fBranchRec("jets"),
 112   fBranchGen(""),
 113   fUseAODInput(kFALSE),
 114   fUseExternalWeightOnly(kFALSE),
 115   fLimitGenJetEta(kFALSE),
 116   fAnalysisType(0),
 117   fExternalWeight(1),
 118   fRecEtaWindow(0.5),
 119   fh1Xsec(0x0),
 120   fh1Trials(0x0),
 121   fh1PtHard(0x0),
 122   fh1PtHardNoW(0x0),
 123   fh1PtHardTrials(0x0),
 124   fh1NGenJets(0x0),
 125   fh1NRecJets(0x0),
 126   fHistList(0x0) ,
 127   ////////////////
 128   fh1JetMultiplicity(0x0) ,
 129   fh2ERecZRec(0x0) ,
 130   fh2EGenZGen(0x0) ,
 131   fh2Efficiency(0x0) ,
 132   fh3EGenERecN(0x0)
 133   ////////////////
 134 {
 135   // Default constructor
 136   for(int ij  = 0;ij<kMaxJets;++ij){
 137     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
 138     fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =0;
 139
 140     fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
 141   }
 142
 143     for(int ic = 0;ic < kMaxCorrelation;ic++){
 144       fhnCorrelation[ic] = 0;
 145     }
 146
 147   DefineOutput(1, TList::Class());
 148 }
 149
 150
 151
 152 Bool_t AliAnalysisTaskJetSpectrum::Notify()
 153 {
 154   //
 155   // Implemented Notify() to read the cross sections
 156   // and number of trials from pyxsec.root
 157   //
 158   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 159   Double_t xsection = 0;
 160   UInt_t   ntrials  = 0;
 161   if(tree){
 162     TFile *curfile = tree->GetCurrentFile();
 163     if (!curfile) {
 164       Error("Notify","No current file");
 165       return kFALSE;
 166     }
 167     if(!fh1Xsec||!fh1Trials){
 168       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 169       return kFALSE;
 170     }
 171
 172     TString fileName(curfile->GetName());
 173     if(fileName.Contains("AliESDs.root")){
 174         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
 175     }
 176     else if(fileName.Contains("AliAOD.root")){
 177         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
 178     }
 179     else if(fileName.Contains("galice.root")){
 180         // for running with galice and kinematics alone...
 181         fileName.ReplaceAll("galice.root", "pyxsec.root");
 182     }
 183     TFile *fxsec = TFile::Open(fileName.Data());
 184     if(!fxsec){
 185       Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
 186       // no a severe condition
 187       return kTRUE;
 188     }
 189     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 190     if(!xtree){
 191       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
 192     }
 193     xtree->SetBranchAddress("xsection",&xsection);
 194     xtree->SetBranchAddress("ntrials",&ntrials);
 195     xtree->GetEntry(0);
 196     fh1Xsec->Fill("<#sigma>",xsection);
 197     fh1Trials->Fill("#sum{ntrials}",ntrials);
 198   }
 199
 200   return kTRUE;
 201 }
 202
 203 void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
 204 {
 205
 206   //
 207   // Create the output container
 208   //
 209
 210
 211   // Connect the AOD
 212
 213   if(fUseAODInput){
 214     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 215     if(!fAOD){
 216       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 217       return;
 218     }
 219     // fethc the header
 220     fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 221     if(!fJetHeaderRec){
 222       Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
 223     }
 224   }
 225   else{
 226     //  assume that the AOD is in the general output...
 227     fAOD  = AODEvent();
 228     if(!fAOD){
 229       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 230       return;
 231     }
 232     fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 233     if(!fJetHeaderRec){
 234       Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
 235     }
 236     else{
 237       if(fDebug>10)fJetHeaderRec->Dump();
 238     }
 239   }
 240
 241
 242
 243   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
 244
 245   OpenFile(1);
 246   if(!fHistList)fHistList = new TList();
 247
 248   Bool_t oldStatus = TH1::AddDirectoryStatus();
 249   TH1::AddDirectory(kFALSE);
 250
 251   //
 252   //  Histogram
 253
 254   const Int_t nBinPt = 100;
 255   Double_t binLimitsPt[nBinPt+1];
 256   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 257     if(iPt == 0){
 258       binLimitsPt[iPt] = 0.0;
 259     }
 260     else {// 1.0
 261       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 262     }
 263   }
 264
 265   const Int_t nBinEta = 26;
 266   Double_t binLimitsEta[nBinEta+1] = {
 267     -1.6,-1.4,-1.2,-1.0,
 268     -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
 269     0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
 270     1.0, 1.2, 1.4, 1.6
 271   };
 272
 273
 274   const Int_t nBinPhi = 30;
 275   Double_t binLimitsPhi[nBinPhi+1];
 276   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 277     if(iPhi==0){
 278       binLimitsPhi[iPhi] = 0;
 279     }
 280     else{
 281       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 282     }
 283   }
 284
 285   const Int_t nBinFrag = 25;
 286
 287
 288   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 289   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 290
 291   fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
 292   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 293
 294   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 295
 296   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 297
 298   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 299
 300   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 301
 302   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 303
 304
 305   for(int ij  = 0;ij<kMaxJets;++ij){
 306     fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
 307     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 308     fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
 309     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 310     fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
 311
 312
 313
 314     fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
 315                              nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 316
 317     fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
 318                              nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
 319
 320     fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
 321                              nBinEta,binLimitsEta,nBinEta,binLimitsEta);
 322
 323
 324     fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}",
 325                                     nBinPt,binLimitsPt,100,-1.0,1.0);
 326     fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
 327                                     nBinPt,binLimitsPt,100,-1.0,1.0);
 328
 329
 330     fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R",
 331                                   nBinPt,binLimitsPt,60,0,6.0);
 332     fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R",
 333                                   nBinPt,binLimitsPt,60,0,6.0);
 334
 335
 336
 337     fh3PtRecGenHard[ij] = new TH3F(Form("fh3PtRecGenHard_j%d",ij), "Pt hard vs. pt gen vs. pt rec;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 338
 339
 340
 341     fh3PtRecGenHardNoW[ij] = new TH3F(Form("fh3PtRecGenHardNoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 342
 343
 344     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
 345                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 346
 347     fh2FragLn[ij] = new TH2F(Form("fh2FragLn_j%d",ij),"Jet Fragmentation Ln;#xi=ln(E_{jet}/E_{i});E_{jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 348                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 349
 350     fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 351                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 352
 353
 354
 355     fh3RecEtaPhiPtNoGen[ij] = new TH3F(Form("fh3RecEtaPhiPtNoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 356                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 357
 358
 359     fh3GenEtaPhiPtNoFound[ij] = new TH3F(Form("fh3GenEtaPhiPtNoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
 360                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 361
 362
 363
 364     fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
 365                                  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 366
 367   }
 368
 369
 370   // tmp histos do not add to the header
 371   TH2F *hCorrPt = new TH2F("fh2PtRecPhiCorrPt","#Delta#phi correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 372   fHistList->Add(hCorrPt);
 373   TH2F *hCorrRanPt = new TH2F("fh2PtRecPhiCorrPtRan","#Delta#phi Random correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 374   fHistList->Add(hCorrRanPt);
 375
 376   TH2F *hCorr = new TH2F("fh2PtRecPhiCorr","#Delta#phi correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 377   fHistList->Add(hCorr);
 378   TH2F *hCorrRan = new TH2F("fh2PtRecPhiCorrRan","#Delta#phi Random correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 379   fHistList->Add(hCorrRan);
 380
 381
 382   /////////////////////////////////////////////////////////////////
 383   fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
 384
 385   fh2ERecZRec   = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
 386   fh2EGenZGen   = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
 387   fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);
 388
 389   fh3EGenERecN  = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100, 0., 250., 100, 0., 250., 51, 0., 50.);
 390
 391   // Response map
 392   //arrays for bin limits
 393   const Int_t nbin[4] = {100, 100, 100, 100};
 394   Double_t vLowEdge[4] = {0.,0.,0.,0.};
 395   Double_t vUpEdge[4] = {250., 250., 1., 1.};
 396
 397   for(int ic = 0;ic < kMaxCorrelation;ic++){
 398     fhnCorrelation[ic]  = new THnSparseF(Form("fhnCorrelation_%d",ic),  "Response Map", 4, nbin, vLowEdge, vUpEdge);
 399     if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
 400     else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
 401   }
 402   const Int_t saveLevel = 3; // large save level more histos
 403   if(saveLevel>0){
 404     fHistList->Add(fh1Xsec);
 405     fHistList->Add(fh1Trials);
 406     fHistList->Add(fh1PtHard);
 407     fHistList->Add(fh1PtHardNoW);
 408     fHistList->Add(fh1PtHardTrials);
 409     fHistList->Add(fh1NGenJets);
 410     fHistList->Add(fh1NRecJets);
 411     ////////////////////////
 412     fHistList->Add(fh1JetMultiplicity);
 413     fHistList->Add(fh2ERecZRec);
 414     fHistList->Add(fh2EGenZGen);
 415     fHistList->Add(fh2Efficiency);
 416     fHistList->Add(fh3EGenERecN);
 417
 418     for(int ic = 0;ic < kMaxCorrelation;++ic){
 419       fHistList->Add(fhnCorrelation[ic]);
 420     }
 421     ////////////////////////
 422     for(int ij  = 0;ij<kMaxJets;++ij){
 423       fHistList->Add(fh1E[ij]);
 424       fHistList->Add(fh1PtRecIn[ij]);
 425       fHistList->Add(fh1PtRecOut[ij]);
 426       fHistList->Add(fh1PtGenIn[ij]);
 427       fHistList->Add(fh1PtGenOut[ij]);
 428       fHistList->Add(fh2PtFGen[ij]);
 429       fHistList->Add(fh2PhiFGen[ij]);
 430       fHistList->Add(fh2EtaFGen[ij]);
 431       fHistList->Add(fh2PtGenDeltaEta[ij]);
 432       fHistList->Add(fh2PtGenDeltaPhi[ij]);
 433       fHistList->Add(fh2PtRecDeltaR[ij]);
 434       fHistList->Add(fh2PtGenDeltaR[ij]);
 435       fHistList->Add(fh3RecEtaPhiPt[ij]);
 436       fHistList->Add(fh3GenEtaPhiPt[ij]);
 437       if(saveLevel>2){
 438         fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
 439         fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
 440       }
 441     }
 442   }
 443
 444   // =========== Switch on Sumw2 for all histos ===========
 445   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 446     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 447     if (h1){
 448       // Printf("%s ",h1->GetName());
 449       h1->Sumw2();
 450       continue;
 451     }
 452     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 453     if(hn)hn->Sumw2();
 454   }
 455
 456   TH1::AddDirectory(oldStatus);
 457
 458 }
 459
 460 void AliAnalysisTaskJetSpectrum::Init()
 461 {
 462   //
 463   // Initialization
 464   //
 465
 466   Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
 467   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
 468
 469 }
 470
 471 void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
 472 {
 473   //
 474   // Execute analysis for current event
 475   //
 476
 477
 478
 479   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 480
 481
 482   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 483
 484   if(!aodH){
 485     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 486     return;
 487   }
 488
 489
 490   //  aodH->SelectEvent(kTRUE);
 491
 492   // ========= These pointers need to be valid in any case =======
 493
 494
 495   /*
 496   AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
 497   if(!jhRec){
 498     Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
 499     return;
 500   }
 501   */
 502   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 503   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 504   if(!aodRecJets){
 505     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 506     return;
 507   }
 508
 509   // ==== General variables needed
 510
 511
 512   // We use statice array, not to fragment the memory
 513   AliAODJet genJets[kMaxJets];
 514   Int_t nGenJets = 0;
 515   AliAODJet recJets[kMaxJets];
 516   Int_t nRecJets = 0;
 517   ///////////////////////////
 518   Int_t nTracks = 0;
 519   //////////////////////////
 520
 521   Double_t eventW = 1;
 522   Double_t ptHard = 0;
 523   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 524
 525   if(fUseExternalWeightOnly){
 526     eventW = fExternalWeight;
 527   }
 528
 529
 530   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 531   if((fAnalysisType&kAnaMC)==kAnaMC){
 532     // this is the part we only use when we have MC information
 533     AliMCEvent* mcEvent = MCEvent();
 534     //    AliStack *pStack = 0;
 535     if(!mcEvent){
 536       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 537       return;
 538     }
 539     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 540     if(!pythiaGenHeader){
 541       return;
 542     }
 543
 544     nTrials = pythiaGenHeader->Trials();
 545     ptHard  = pythiaGenHeader->GetPtHard();
 546     int iProcessType = pythiaGenHeader->ProcessType();
 547     // 11 f+f -> f+f
 548     // 12 f+barf -> f+barf
 549     // 13 f+barf -> g+g
 550     // 28 f+g -> f+g
 551     // 53 g+g -> f+barf
 552     // 68 g+g -> g+g
 553     /*
 554     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 555     //    if(iProcessType != 13 && iProcessType != 68){ // allow only glue
 556     if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
 557     //    if(iProcessType != 28){ // allow only -> f+g
 558       PostData(1, fHistList);
 559       return;
 560     }
 561     */
 562     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 563
 564     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 565
 566     if(!fUseExternalWeightOnly){
 567         // case were we combine more than one p_T hard bin...
 568     }
 569
 570     // fetch the pythia generated jets only to be used here
 571     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 572     AliAODJet pythiaGenJets[kMaxJets];
 573     Int_t iCount = 0;
 574     for(int ip = 0;ip < nPythiaGenJets;++ip){
 575       if(iCount>=kMaxJets)continue;
 576       Float_t p[4];
 577       pythiaGenHeader->TriggerJet(ip,p);
 578       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 579
 580       if(fLimitGenJetEta){
 581         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 582            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 583       }
 584
 585
 586       if(fBranchGen.Length()==0){
 587         // if we have MC particles and we do not read from the aod branch
 588         // use the pythia jets
 589         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 590       }
 591       iCount++;
 592     }
 593     if(fBranchGen.Length()==0)nGenJets = iCount;
 594
 595   }// (fAnalysisType&kMC)==kMC)
 596
 597   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 598   fh1PtHard->Fill(ptHard,eventW);
 599   fh1PtHardNoW->Fill(ptHard,1);
 600   fh1PtHardTrials->Fill(ptHard,nTrials);
 601
 602   // If we set a second branch for the input jets fetch this
 603   if(fBranchGen.Length()>0){
 604     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 605     if(aodGenJets){
 606       Int_t iCount = 0;
 607       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 608         if(iCount>=kMaxJets)continue;
 609         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 610         if(!tmp)continue;
 611         if(fLimitGenJetEta){
 612           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 613              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 614         }
 615         genJets[iCount] = *tmp;
 616         iCount++;
 617       }
 618       nGenJets = iCount;
 619     }
 620     else{
 621       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 622     }
 623   }
 624
 625   fh1NGenJets->Fill(nGenJets);
 626   // We do not want to exceed the maximum jet number
 627   nGenJets = TMath::Min(nGenJets,kMaxJets);
 628
 629   // Fetch the reconstructed jets...
 630
 631
 632   nRecJets = aodRecJets->GetEntries();
 633   fh1NRecJets->Fill(nRecJets);
 634   nRecJets = TMath::Min(nRecJets,kMaxJets);
 635   //////////////////////////////////////////
 636   nTracks  = fAOD->GetNumberOfTracks();
 637   ///////////////////////////////////////////
 638
 639   for(int ir = 0;ir < nRecJets;++ir){
 640     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 641     if(!tmp)continue;
 642     recJets[ir] = *tmp;
 643   }
 644
 645
 646   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 647   // Relate the jets
 648   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 649   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 650
 651   for(int i = 0;i<kMaxJets;++i){
 652     iGenIndex[i] = iRecIndex[i] = -1;
 653   }
 654
 655
 656   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 657                  iGenIndex,iRecIndex,fDebug);
 658   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 659
 660   if(fDebug){
 661     for(int i = 0;i<kMaxJets;++i){
 662       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 663       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 664     }
 665   }
 666
 667   // loop over reconstructed jets
 668   for(int ir = 0;ir < nRecJets;++ir){
 669     Double_t ptRec = recJets[ir].Pt();
 670     Double_t phiRec = recJets[ir].Phi();
 671     Double_t phiRecRan = TMath::Pi()*gRandom->Rndm(); // better take real jet axis from previous events (TPC acceptance in phi)
 672     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 673     Double_t etaRec = recJets[ir].Eta();
 674     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 675     fh1E[ir]->Fill(recJets[ir].E(),eventW);
 676     fh1PtRecIn[ir]->Fill(ptRec,eventW);
 677     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 678     for(int irr = ir+1;irr<nRecJets;irr++){
 679       fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
 680     }
 681     // Fill Correlation
 682     Int_t ig = iGenIndex[ir];
 683     if(ig>=0 && ig<nGenJets){
 684       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 685       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 686       fh1PtRecOut[ir]->Fill(ptRec,eventW);
 687       Double_t ptGen  = genJets[ig].Pt();
 688       Double_t phiGen = genJets[ig].Phi();
 689       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 690       Double_t etaGen = genJets[ig].Eta();
 691
 692       //
 693       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
 694       //
 695
 696       if(TMath::Abs(etaRec)<fRecEtaWindow){
 697
 698       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
 699       fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
 700       fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
 701       fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
 702       fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
 703       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
 704       fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
 705       /////////////////////////////////////////////////////
 706
 707       //      Double_t eRec = recJets[ir].E();
 708       //      Double_t eGen = genJets[ig].E();
 709       // CKB use p_T not Energy
 710       // TODO recname variabeles and histos
 711       Double_t eRec = recJets[ir].E();
 712       Double_t eGen = genJets[ig].E();
 713
 714       fh2Efficiency->Fill(eGen, eRec/eGen);
 715
 716       if (eGen>=0. && eGen<=250.){
 717         Double_t eLeading = -1;
 718         Double_t ptleading = -1;
 719         Int_t nPart=0;
 720         // loop over tracks
 721         for (Int_t it = 0; it< nTracks; it++){
 722           //      if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
 723           // find leading particle
 724           //  if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
 725           // TODO implement esd filter flag to be the same as in the jet finder
 726           // allow also for MC particles...
 727           Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
 728           if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
 729             eLeading  = fAOD->GetTrack(it)->E();
 730             ptleading = fAOD->GetTrack(it)->Pt();
 731           }
 732           //          if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties
 733           if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
 734             nPart++;
 735
 736
 737           // correlate jet axis of leading jet with particles
 738           if(ir==0){
 739             Float_t phi =  fAOD->GetTrack(it)->Phi();
 740             Float_t dPhi = phi - phiRec;
 741             if(dPhi>TMath::Pi()/1.5)dPhi = dPhi - 2.*TMath::Pi();
 742             if(dPhi<(-0.5*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 743             Float_t dPhiRan = phi - phiRecRan;
 744             if(dPhiRan>TMath::Pi()/1.5)dPhiRan = dPhiRan - 2.*TMath::Pi();
 745             if(dPhiRan<(-0.5*TMath::Pi()))dPhiRan = dPhiRan + 2.*TMath::Pi();
 746             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorr"))->Fill(ptRec,dPhi);
 747             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrRan"))->Fill(ptRec,dPhiRan);
 748             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPt"))->Fill(ptRec,dPhi,fAOD->GetTrack(it)->Pt());
 749             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPtRan"))->Fill(ptRec,dPhiRan,fAOD->GetTrack(it)->Pt());
 750
 751           }
 752         }
 753         if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 754
 755         // fill Response Map (4D histogram) and Energy vs z distributions
 756         Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
 757         fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
 758         fh2EGenZGen->Fill(var[0],var[2]);
 759         fh1JetMultiplicity->Fill(nPart);
 760         fh3EGenERecN->Fill(eGen, eRec, nPart);
 761         for(int ic = 0;ic <kMaxCorrelation;ic++){
 762           if (nPart<=fgkJetNpartCut[ic]){ // is this corrected for CKB
 763             fhnCorrelation[ic]->Fill(var);
 764             break;
 765           }
 766         }
 767       }
 768
 769     }// if etarec in window
 770
 771     }
 772     ////////////////////////////////////////////////////
 773     else{
 774       fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 775     }
 776   }// loop over reconstructed jets
 777
 778
 779   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 780   for(int ig = 0;ig < nGenJets;++ig){
 781     Double_t ptGen = genJets[ig].Pt();
 782     // Fill Correlation
 783     Double_t phiGen = genJets[ig].Phi();
 784     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 785     Double_t etaGen = genJets[ig].Eta();
 786     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 787     fh1PtGenIn[ig]->Fill(ptGen,eventW);
 788     for(int igg = ig+1;igg<nGenJets;igg++){
 789       fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
 790     }
 791     Int_t ir = iRecIndex[ig];
 792     if(ir>=0&&ir<nRecJets){
 793       fh1PtGenOut[ig]->Fill(ptGen,eventW);
 794     }
 795     else{
 796       fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 797     }
 798   }// loop over reconstructed jets
 799
 800   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 801   PostData(1, fHistList);
 802 }
 803
 804 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
 805 {
 806 // Terminate analysis
 807 //
 808     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
 809 }