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