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 <TSystem.h>
  25 #include <TInterpreter.h>
  26 #include <TChain.h>
  27 #include <TFile.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TH3F.h>
  31 #include <TProfile.h>
  32 #include <TList.h>
  33 #include <TLorentzVector.h>
  34 #include <TClonesArray.h>
  35 #include  "TDatabasePDG.h"
  36
  37 #include "AliAnalysisTaskJetSpectrum.h"
  38 #include "AliAnalysisManager.h"
  39 #include "AliJetFinder.h"
  40 #include "AliJetHeader.h"
  41 #include "AliJetReader.h"
  42 #include "AliJetReaderHeader.h"
  43 #include "AliUA1JetHeaderV1.h"
  44 #include "AliJet.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliAODHandler.h"
  48 #include "AliAODTrack.h"
  49 #include "AliAODJet.h"
  50 #include "AliMCEventHandler.h"
  51 #include "AliMCEvent.h"
  52 #include "AliStack.h"
  53 #include "AliGenPythiaEventHeader.h"
  54 #include "AliJetKineReaderHeader.h"
  55 #include "AliGenCocktailEventHeader.h"
  56 #include "AliInputEventHandler.h"
  57
  58 #include "AliAnalysisHelperJetTasks.h"
  59
  60 ClassImp(AliAnalysisTaskJetSpectrum)
  61
  62 const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
  63
  64 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
  65   fJetHeaderRec(0x0),
  66   fJetHeaderGen(0x0),
  67   fAOD(0x0),
  68   fBranchRec("jets"),
  69   fConfigRec("ConfigJets.C"),
  70   fBranchGen(""),
  71   fConfigGen(""),
  72   fUseAODInput(kFALSE),
  73   fUseExternalWeightOnly(kFALSE),
  74   fLimitGenJetEta(kFALSE),
  75   fAnalysisType(0),
  76   fExternalWeight(1),
  77   fh1Xsec(0x0),
  78   fh1Trials(0x0),
  79   fh1PtHard(0x0),
  80   fh1PtHardNoW(0x0),
  81   fh1PtHardTrials(0x0),
  82   fh1NGenJets(0x0),
  83   fh1NRecJets(0x0),
  84   fHistList(0x0) ,
  85   ////////////////
  86   fh1JetMultiplicity(0x0) ,
  87   fh2ERecZRec(0x0) ,
  88   fh2EGenZGen(0x0) ,
  89   fh2Efficiency(0x0) ,
  90   fh3EGenERecN(0x0)
  91   ////////////////
  92 {
  93   // Default constructor
  94     for(int ij  = 0;ij<kMaxJets;++ij){
  95       fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
  96       fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] =  fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
  97       fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
  98     }
  99     for(int ic = 0;ic < kMaxCorrelation;ic++){
 100       fhnCorrelation[ic] = 0;
 101     }
 102
 103 }
 104
 105 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
 106   AliAnalysisTaskSE(name),
 107   fJetHeaderRec(0x0),
 108   fJetHeaderGen(0x0),
 109   fAOD(0x0),
 110   fBranchRec("jets"),
 111   fConfigRec("ConfigJets.C"),
 112   fBranchGen(""),
 113   fConfigGen(""),
 114   fUseAODInput(kFALSE),
 115   fUseExternalWeightOnly(kFALSE),
 116   fLimitGenJetEta(kFALSE),
 117   fAnalysisType(0),
 118   fExternalWeight(1),
 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] = 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;
 262     }
 263   }
 264   const Int_t nBinEta = 22;
 265   Double_t binLimitsEta[nBinEta+1];
 266   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 267     if(iEta==0){
 268       binLimitsEta[iEta] = -2.2;
 269     }
 270     else{
 271       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.2;
 272     }
 273   }
 274
 275   const Int_t nBinPhi = 30;
 276   Double_t binLimitsPhi[nBinPhi+1];
 277   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 278     if(iPhi==0){
 279       binLimitsPhi[iPhi] = 0;
 280     }
 281     else{
 282       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 283     }
 284   }
 285
 286   const Int_t nBinFrag = 25;
 287
 288
 289   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 290   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 291
 292   fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
 293   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 294
 295   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 296
 297   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 298
 299   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 300
 301   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 302
 303   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 304
 305
 306   for(int ij  = 0;ij<kMaxJets;++ij){
 307     fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
 308     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 309     fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
 310     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 311     fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
 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
 331
 332
 333     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);
 334
 335
 336
 337     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);
 338
 339
 340     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
 341                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 342
 343     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",
 344                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 345
 346     fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 347                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 348
 349
 350
 351     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)",
 352                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 353
 354
 355     fh3GenEtaPhiPtNoFound[ij] = new TH3F(Form("fh3GenEtaPhiPtNoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
 356                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 357
 358
 359
 360     fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 361                                  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 362
 363   }
 364
 365   /////////////////////////////////////////////////////////////////
 366   fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
 367
 368   fh2ERecZRec   = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
 369   fh2EGenZGen   = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
 370   fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);
 371
 372   fh3EGenERecN  = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100, 0., 250., 100, 0., 250., 51, 0., 50.);
 373
 374   // Response map
 375   //arrays for bin limits
 376   const Int_t nbin[4] = {100, 100, 100, 100};
 377   Double_t vLowEdge[4] = {0.,0.,0.,0.};
 378   Double_t vUpEdge[4] = {250., 250., 1., 1.};
 379
 380   for(int ic = 0;ic < kMaxCorrelation;ic++){
 381     fhnCorrelation[ic]  = new THnSparseF(Form("fhnCorrelation_%d",ic),  "Response Map", 4, nbin, vLowEdge, vUpEdge);
 382     if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
 383     else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
 384   }
 385   const Int_t saveLevel = 3; // large save level more histos
 386   if(saveLevel>0){
 387     fHistList->Add(fh1Xsec);
 388     fHistList->Add(fh1Trials);
 389     fHistList->Add(fh1PtHard);
 390     fHistList->Add(fh1PtHardNoW);
 391     fHistList->Add(fh1PtHardTrials);
 392     fHistList->Add(fh1NGenJets);
 393     fHistList->Add(fh1NRecJets);
 394     ////////////////////////
 395     fHistList->Add(fh1JetMultiplicity);
 396     fHistList->Add(fh2ERecZRec);
 397     fHistList->Add(fh2EGenZGen);
 398     fHistList->Add(fh2Efficiency);
 399     fHistList->Add(fh3EGenERecN);
 400
 401     for(int ic = 0;ic < kMaxCorrelation;++ic){
 402       fHistList->Add(fhnCorrelation[ic]);
 403     }
 404     ////////////////////////
 405     for(int ij  = 0;ij<kMaxJets;++ij){
 406       fHistList->Add(fh1E[ij]);
 407       fHistList->Add(fh1PtRecIn[ij]);
 408       fHistList->Add(fh1PtRecOut[ij]);
 409       fHistList->Add(fh1PtGenIn[ij]);
 410       fHistList->Add(fh1PtGenOut[ij]);
 411       fHistList->Add(fh2PtFGen[ij]);
 412       fHistList->Add(fh2PhiFGen[ij]);
 413       fHistList->Add(fh2EtaFGen[ij]);
 414       fHistList->Add(fh2PtGenDeltaEta[ij]);
 415       fHistList->Add(fh2PtGenDeltaPhi[ij]);
 416       fHistList->Add(fh3RecEtaPhiPt[ij]);
 417       fHistList->Add(fh3GenEtaPhiPt[ij]);
 418       if(saveLevel>2){
 419         fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
 420         fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
 421       }
 422     }
 423   }
 424
 425   // =========== Switch on Sumw2 for all histos ===========
 426   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 427     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 428     if (h1){
 429       // Printf("%s ",h1->GetName());
 430       h1->Sumw2();
 431       continue;
 432     }
 433     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 434     if(hn)hn->Sumw2();
 435   }
 436
 437   TH1::AddDirectory(oldStatus);
 438
 439 }
 440
 441 void AliAnalysisTaskJetSpectrum::Init()
 442 {
 443   //
 444   // Initialization
 445   //
 446
 447   Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
 448   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
 449
 450 }
 451
 452 void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
 453 {
 454   //
 455   // Execute analysis for current event
 456   //
 457
 458
 459
 460   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 461
 462
 463   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 464
 465   if(!aodH){
 466     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 467     return;
 468   }
 469
 470
 471   //  aodH->SelectEvent(kTRUE);
 472
 473   // ========= These pointers need to be valid in any case =======
 474
 475
 476   /*
 477   AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
 478   if(!jhRec){
 479     Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
 480     return;
 481   }
 482   */
 483   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 484   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 485   if(!aodRecJets){
 486     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 487     return;
 488   }
 489
 490   // ==== General variables needed
 491
 492
 493   // We use statice array, not to fragment the memory
 494   AliAODJet genJets[kMaxJets];
 495   Int_t nGenJets = 0;
 496   AliAODJet recJets[kMaxJets];
 497   Int_t nRecJets = 0;
 498   ///////////////////////////
 499   Int_t nTracks = 0;
 500   //////////////////////////
 501
 502   Double_t eventW = 1;
 503   Double_t ptHard = 0;
 504   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 505
 506   if(fUseExternalWeightOnly){
 507     eventW = fExternalWeight;
 508   }
 509
 510
 511   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 512   if((fAnalysisType&kAnaMC)==kAnaMC){
 513     // this is the part we only use when we have MC information
 514     AliMCEvent* mcEvent = MCEvent();
 515     //    AliStack *pStack = 0;
 516     if(!mcEvent){
 517       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 518       return;
 519     }
 520     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 521     if(!pythiaGenHeader){
 522       return;
 523     }
 524
 525     nTrials = pythiaGenHeader->Trials();
 526     ptHard  = pythiaGenHeader->GetPtHard();
 527
 528
 529     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 530
 531     if(!fUseExternalWeightOnly){
 532         // case were we combine more than one p_T hard bin...
 533     }
 534
 535     // fetch the pythia generated jets only to be used here
 536     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 537     AliAODJet pythiaGenJets[kMaxJets];
 538     Int_t iCount = 0;
 539     for(int ip = 0;ip < nPythiaGenJets;++ip){
 540       if(iCount>=kMaxJets)continue;
 541       Float_t p[4];
 542       pythiaGenHeader->TriggerJet(ip,p);
 543       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 544
 545       if(fLimitGenJetEta){
 546         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 547            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 548       }
 549
 550
 551       if(fBranchGen.Length()==0){
 552         // if we have MC particles and we do not read from the aod branch
 553         // use the pythia jets
 554         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 555       }
 556       iCount++;
 557     }
 558     if(fBranchGen.Length()==0)nGenJets = iCount;
 559
 560   }// (fAnalysisType&kMC)==kMC)
 561
 562   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 563   fh1PtHard->Fill(ptHard,eventW);
 564   fh1PtHardNoW->Fill(ptHard,1);
 565   fh1PtHardTrials->Fill(ptHard,nTrials);
 566
 567   // If we set a second branch for the input jets fetch this
 568   if(fBranchGen.Length()>0){
 569     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 570     if(aodGenJets){
 571       Int_t iCount = 0;
 572       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 573         if(iCount>=kMaxJets)continue;
 574         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 575         if(!tmp)continue;
 576         if(fLimitGenJetEta){
 577           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 578              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 579         }
 580         genJets[iCount] = *tmp;
 581         iCount++;
 582       }
 583       nGenJets = iCount;
 584     }
 585     else{
 586       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 587     }
 588   }
 589
 590   fh1NGenJets->Fill(nGenJets);
 591   // We do not want to exceed the maximum jet number
 592   nGenJets = TMath::Min(nGenJets,kMaxJets);
 593
 594   // Fetch the reconstructed jets...
 595
 596
 597   nRecJets = aodRecJets->GetEntries();
 598   fh1NRecJets->Fill(nRecJets);
 599   nRecJets = TMath::Min(nRecJets,kMaxJets);
 600   //////////////////////////////////////////
 601   nTracks  = fAOD->GetNumberOfTracks();
 602   ///////////////////////////////////////////
 603
 604   for(int ir = 0;ir < nRecJets;++ir){
 605     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 606     if(!tmp)continue;
 607     recJets[ir] = *tmp;
 608   }
 609
 610
 611   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 612   // Relate the jets
 613   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 614   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 615
 616   for(int i = 0;i<kMaxJets;++i){
 617     iGenIndex[i] = iRecIndex[i] = -1;
 618   }
 619
 620
 621   GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 622                  iGenIndex,iRecIndex,fDebug);
 623   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 624
 625   if(fDebug){
 626     for(int i = 0;i<kMaxJets;++i){
 627       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 628       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 629     }
 630   }
 631
 632   // loop over reconstructed jets
 633   for(int ir = 0;ir < nRecJets;++ir){
 634     Double_t ptRec = recJets[ir].Pt();
 635     Double_t phiRec = recJets[ir].Phi();
 636     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 637     Double_t etaRec = recJets[ir].Eta();
 638     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 639     fh1E[ir]->Fill(recJets[ir].E(),eventW);
 640     fh1PtRecIn[ir]->Fill(ptRec,eventW);
 641     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 642     // Fill Correlation
 643     Int_t ig = iGenIndex[ir];
 644     if(ig>=0 && ig<nGenJets){
 645       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 646       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 647       fh1PtRecOut[ir]->Fill(ptRec,eventW);
 648       Double_t ptGen  = genJets[ig].Pt();
 649       Double_t phiGen = genJets[ig].Phi();
 650       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 651       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 652       Double_t etaGen = genJets[ig].Eta();
 653       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 654       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
 655       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 656       fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
 657       fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
 658       fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
 659       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 660       fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
 661       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 662       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
 663       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 664       fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
 665       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 666       /////////////////////////////////////////////////////
 667       Double_t eRec = recJets[ir].E();
 668       Double_t eGen = genJets[ig].E();
 669       fh2Efficiency->Fill(eGen, eRec/eGen);
 670       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 671
 672       if (eGen>=0. && eGen<=250.){
 673         Double_t eLeading = -1;
 674         Double_t ptleading = -1;
 675         Int_t nPart=0;
 676         // loop over tracks
 677         for (Int_t it = 0; it< nTracks; it++){
 678           if (fAOD->GetTrack(it)->E() > eGen) continue;
 679           Double_t phiTrack = fAOD->GetTrack(it)->Phi();
 680           if (phiTrack<0) phiTrack+=TMath::Pi()*2.;
 681           Double_t etaTrack = fAOD->GetTrack(it)->Eta();
 682           Float_t deta  = etaRec - etaTrack;
 683           Float_t dphi  = TMath::Abs(phiRec - phiTrack);
 684           Float_t r = TMath::Sqrt(deta*deta + dphi*dphi);
 685           // find leading particle
 686           if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
 687             eLeading  = fAOD->GetTrack(it)->E();
 688             ptleading = fAOD->GetTrack(it)->Pt();
 689           }
 690           if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7)
 691             nPart++;
 692         }
 693         if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 694
 695         // fill Response Map (4D histogram) and Energy vs z distributions
 696         Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
 697         fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
 698         fh2EGenZGen->Fill(var[0],var[2]);
 699         fh1JetMultiplicity->Fill(nPart);
 700         fh3EGenERecN->Fill(eGen, eRec, nPart);
 701         for(int ic = 0;ic <kMaxCorrelation;ic++){
 702           if (nPart<=fgkJetNpartCut[ic]){
 703             fhnCorrelation[ic]->Fill(var);
 704             break;
 705           }
 706         }
 707       }
 708     }
 709     ////////////////////////////////////////////////////
 710     else{
 711       fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 712     }
 713   }// loop over reconstructed jets
 714
 715
 716   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 717   for(int ig = 0;ig < nGenJets;++ig){
 718     Double_t ptGen = genJets[ig].Pt();
 719     // Fill Correlation
 720     Double_t phiGen = genJets[ig].Phi();
 721     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 722     Double_t etaGen = genJets[ig].Eta();
 723     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 724     fh1PtGenIn[ig]->Fill(ptGen,eventW);
 725     Int_t ir = iRecIndex[ig];
 726     if(ir>=0&&ir<nRecJets){
 727       fh1PtGenOut[ig]->Fill(ptGen,eventW);
 728     }
 729     else{
 730       fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 731     }
 732   }// loop over reconstructed jets
 733
 734   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 735   PostData(1, fHistList);
 736 }
 737
 738 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
 739 {
 740 // Terminate analysis
 741 //
 742     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
 743 }
 744
 745
 746 void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
 747                                                 AliAODJet *recJets,const Int_t &nRecJets,
 748                                                 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
 749
 750   //
 751   // Relate the two input jet Arrays
 752   //
 753
 754   //
 755   // The association has to be unique
 756   // So check in two directions
 757   // find the closest rec to a gen
 758   // and check if there is no other rec which is closer
 759   // Caveat: Close low energy/split jets may disturb this correlation
 760
 761   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
 762   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
 763   // in the end we have something like this
 764   //    1   2   3
 765   // ------------
 766   // a| 3   2   0
 767   // b| 0   1   0
 768   // c| 0   0   3
 769   // d| 0   0   1
 770   // e| 0   0   1
 771   // Topology
 772   //   1     2
 773   //     a         b
 774   //
 775   //  d      c
 776   //        3     e
 777   // Only entries with "3" match from both sides
 778
 779   const int kMode = 3;
 780
 781   Int_t iFlag[kMaxJets][kMaxJets];
 782
 783
 784
 785   for(int i = 0;i < kMaxJets;++i){
 786     iRecIndex[i] = -1;
 787     iGenIndex[i] = -1;
 788     for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
 789   }
 790
 791   if(nRecJets==0)return;
 792   if(nGenJets==0)return;
 793
 794   const Float_t maxDist = 1.0;
 795   // find the closest distance to the generated
 796   for(int ig = 0;ig<nGenJets;++ig){
 797     Float_t dist = maxDist;
 798     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());
 799     for(int ir = 0;ir<nRecJets;++ir){
 800       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
 801       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());
 802       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
 803       if(dR<dist){
 804         iRecIndex[ig] = ir;
 805         dist = dR;
 806       }
 807     }
 808     if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
 809     // reset...
 810     iRecIndex[ig] = -1;
 811   }
 812   // other way around
 813   for(int ir = 0;ir<nRecJets;++ir){
 814     Float_t dist = maxDist;
 815     for(int ig = 0;ig<nGenJets;++ig){
 816       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
 817       if(dR<dist){
 818         iGenIndex[ir] = ig;
 819         dist = dR;
 820       }
 821     }
 822     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
 823     // reset...
 824     iGenIndex[ir] = -1;
 825   }
 826
 827   // check for "true" correlations
 828
 829   if(iDebug>1)Printf(">>>>>> Matrix");
 830
 831   for(int ig = 0;ig<nGenJets;++ig){
 832     for(int ir = 0;ir<nRecJets;++ir){
 833       // Print
 834       if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
 835
 836       if(kMode==3){
 837         // we have a uniqie correlation
 838         if(iFlag[ig][ir]==3){
 839           iGenIndex[ir] = ig;
 840           iRecIndex[ig] = ir;
 841         }
 842       }
 843       else{
 844         // we just take the correlation from on side
 845         if((iFlag[ig][ir]&2)==2){
 846           iGenIndex[ir] = ig;
 847         }
 848         if((iFlag[ig][ir]&1)==1){
 849           iRecIndex[ig] = ir;
 850         }
 851       }
 852     }
 853     if(iDebug>1)printf("\n");
 854   }
 855 }
 856
 857