PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx

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