PWG4/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 <TList.h>
  27 #include <TLorentzVector.h>
  28 #include <TClonesArray.h>
  29 #include  "TDatabasePDG.h"
  30
  31 #include "AliAnalysisTaskJetSpectrum.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliJetFinder.h"
  34 #include "AliJetHeader.h"
  35 #include "AliJetReader.h"
  36 #include "AliJetReaderHeader.h"
  37 #include "AliUA1JetHeaderV1.h"
  38 #include "AliJet.h"
  39 #include "AliESDEvent.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODHandler.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODJet.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliStack.h"
  47 #include "AliGenPythiaEventHeader.h"
  48 #include "AliJetKineReaderHeader.h"
  49 #include "AliGenCocktailEventHeader.h"
  50 #include "AliInputEventHandler.h"
  51
  52 #include "AliAnalysisHelperJetTasks.h"
  53
  54 ClassImp(AliAnalysisTaskJetSpectrum)
  55
  56 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
  57   fJetHeaderRec(0x0),
  58   fJetHeaderGen(0x0),
  59   fAOD(0x0),
  60   fBranchRec("jets"),
  61   fConfigRec("ConfigJets.C"),
  62   fBranchGen(""),
  63   fConfigGen(""),
  64   fUseAODInput(kFALSE),
  65   fUseExternalWeightOnly(kFALSE),
  66   fLimitGenJetEta(kFALSE),
  67   fAnalysisType(0),
  68   fExternalWeight(1),
  69   fh1PtHard(0x0),
  70   fh1PtHard_NoW(0x0),
  71   fh1PtHard_Trials(0x0),
  72   fh1NGenJets(0x0),
  73   fh1NRecJets(0x0),
  74   fHistList(0x0)
  75 {
  76   // Default constructor
  77     for(int ij  = 0;ij<kMaxJets;++ij){
  78     fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
  79     fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
  80     fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
  81   }
  82
  83 }
  84
  85 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
  86   AliAnalysisTaskSE(name),
  87   fJetHeaderRec(0x0),
  88   fJetHeaderGen(0x0),
  89   fAOD(0x0),
  90   fBranchRec("jets"),
  91   fConfigRec("ConfigJets.C"),
  92   fBranchGen(""),
  93   fConfigGen(""),
  94   fUseAODInput(kFALSE),
  95   fUseExternalWeightOnly(kFALSE),
  96   fLimitGenJetEta(kFALSE),
  97   fAnalysisType(0),
  98   fExternalWeight(1),
  99   fh1PtHard(0x0),
 100   fh1PtHard_NoW(0x0),
 101   fh1PtHard_Trials(0x0),
 102   fh1NGenJets(0x0),
 103   fh1NRecJets(0x0),
 104   fHistList(0x0)
 105 {
 106   // Default constructor
 107   for(int ij  = 0;ij<kMaxJets;++ij){
 108     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
 109     fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
 110
 111     fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
 112   }
 113
 114   DefineOutput(1, TList::Class());
 115 }
 116
 117
 118
 119
 120 void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
 121 {
 122
 123   //
 124   // Create the output container
 125   //
 126
 127
 128   // Connect the AOD
 129
 130   if(fUseAODInput){
 131     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 132     if(!fAOD){
 133       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 134       return;
 135     }
 136     // fethc the header
 137     fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 138     if(!fJetHeaderRec){
 139       Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
 140     }
 141   }
 142   else{
 143     //  assume that the AOD is in the general output...
 144     fAOD  = AODEvent();
 145     if(!fAOD){
 146       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 147       return;
 148     }
 149     fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 150     if(!fJetHeaderRec){
 151       Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
 152     }
 153     else{
 154       if(fDebug>10)fJetHeaderRec->Dump();
 155     }
 156   }
 157
 158
 159
 160   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
 161
 162   OpenFile(1);
 163   if(!fHistList)fHistList = new TList();
 164
 165   Bool_t oldStatus = TH1::AddDirectoryStatus();
 166   TH1::AddDirectory(kFALSE);
 167
 168   //
 169   //  Histogram
 170
 171   const Int_t nBinPt = 100;
 172   Double_t binLimitsPt[nBinPt+1];
 173   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 174     if(iPt == 0){
 175       binLimitsPt[iPt] = 0.0;
 176     }
 177     else {// 1.0
 178       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2;
 179     }
 180   }
 181   const Int_t nBinEta = 22;
 182   Double_t binLimitsEta[nBinEta+1];
 183   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 184     if(iEta==0){
 185       binLimitsEta[iEta] = -2.2;
 186     }
 187     else{
 188       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.2;
 189     }
 190   }
 191
 192   const Int_t nBinPhi = 90;
 193   Double_t binLimitsPhi[nBinPhi+1];
 194   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 195     if(iPhi==0){
 196       binLimitsPhi[iPhi] = 0;
 197     }
 198     else{
 199       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 200     }
 201   }
 202
 203   const Int_t nBinFrag = 25;
 204
 205
 206   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 207
 208   fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 209
 210   fh1PtHard_Trials = new TH1F("fh1PtHard_Trials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 211
 212   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 213
 214   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 215
 216
 217   for(int ij  = 0;ij<kMaxJets;++ij){
 218     fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
 219     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 220     fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
 221     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 222     fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
 223
 224
 225     fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
 226                              nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 227
 228
 229
 230
 231
 232
 233     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);
 234
 235
 236
 237     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);
 238
 239
 240     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
 241                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 242
 243     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",
 244                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 245
 246     fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 247                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 248
 249
 250
 251     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)",
 252                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 253
 254
 255     fh3GenEtaPhiPt_NoFound[ij] = new TH3F(Form("fh3GenEtaPhiPt_NoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 256                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 257
 258
 259
 260     fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 261                                  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 262
 263   }
 264
 265   const Int_t saveLevel = 2; // large save level more histos
 266
 267   if(saveLevel>0){
 268     fHistList->Add(fh1PtHard);
 269     fHistList->Add(fh1PtHard_NoW);
 270     fHistList->Add(fh1PtHard_Trials);
 271     fHistList->Add(fh1NGenJets);
 272     fHistList->Add(fh1NRecJets);
 273     for(int ij  = 0;ij<kMaxJets;++ij){
 274       fHistList->Add(fh1E[ij]);
 275       fHistList->Add(fh1PtRecIn[ij]);
 276       fHistList->Add(fh1PtRecOut[ij]);
 277       fHistList->Add(fh1PtGenIn[ij]);
 278       fHistList->Add(fh1PtGenOut[ij]);
 279       fHistList->Add(fh2PtFGen[ij]);
 280       if(saveLevel>2){
 281         fHistList->Add(fh3RecEtaPhiPt[ij]);
 282         fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
 283         fHistList->Add(fh3GenEtaPhiPt_NoFound[ij]);
 284         fHistList->Add(fh3GenEtaPhiPt[ij]);
 285       }
 286     }
 287   }
 288
 289   // =========== Switch on Sumw2 for all histos ===========
 290   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 291     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 292     if (h1){
 293       // Printf("%s ",h1->GetName());
 294       h1->Sumw2();
 295     }
 296   }
 297
 298   TH1::AddDirectory(oldStatus);
 299
 300 }
 301
 302 void AliAnalysisTaskJetSpectrum::Init()
 303 {
 304   //
 305   // Initialization
 306   //
 307
 308   Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
 309   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
 310
 311 }
 312
 313 void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
 314 {
 315   //
 316   // Execute analysis for current event
 317   //
 318
 319
 320
 321   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 322
 323
 324   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 325
 326   if(!aodH){
 327     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 328     return;
 329   }
 330
 331
 332   //  aodH->SelectEvent(kTRUE);
 333
 334   // ========= These pointers need to be valid in any case =======
 335
 336
 337   /*
 338   AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
 339   if(!jhRec){
 340     Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
 341     return;
 342   }
 343   */
 344   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 345   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 346   if(!aodRecJets){
 347     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 348     return;
 349   }
 350
 351   // ==== General variables needed
 352
 353
 354   // We use statice array, not to fragment the memory
 355   AliAODJet genJets[kMaxJets];
 356   Int_t nGenJets = 0;
 357   AliAODJet recJets[kMaxJets];
 358   Int_t nRecJets = 0;
 359
 360   Double_t eventW = 1;
 361   Double_t ptHard = 0;
 362   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 363
 364   if(fUseExternalWeightOnly){
 365     eventW = fExternalWeight;
 366   }
 367
 368
 369   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 370   if((fAnalysisType&kAnaMC)==kAnaMC){
 371     // this is the part we only use when we have MC information
 372     AliMCEvent* mcEvent = MCEvent();
 373     //    AliStack *pStack = 0;
 374     if(!mcEvent){
 375       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 376       return;
 377     }
 378     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 379     if(!pythiaGenHeader){
 380       return;
 381     }
 382
 383     nTrials = pythiaGenHeader->Trials();
 384     ptHard  = pythiaGenHeader->GetPtHard();
 385
 386
 387     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 388
 389     if(!fUseExternalWeightOnly){
 390         // case were we combine more than one p_T hard bin...
 391     }
 392
 393     // fetch the pythia generated jets only to be used here
 394     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 395     AliAODJet pythiaGenJets[kMaxJets];
 396     Int_t iCount = 0;
 397     for(int ip = 0;ip < nPythiaGenJets;++ip){
 398       if(iCount>=kMaxJets)continue;
 399       Float_t p[4];
 400       pythiaGenHeader->TriggerJet(ip,p);
 401       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 402
 403       if(fLimitGenJetEta){
 404         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 405            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 406       }
 407
 408
 409       if(fBranchGen.Length()==0){
 410         // if we have MC particles and we do not read from the aod branch
 411         // use the pythia jets
 412         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 413       }
 414       iCount++;
 415     }
 416     if(fBranchGen.Length()==0)nGenJets = iCount;
 417
 418   }// (fAnalysisType&kMC)==kMC)
 419
 420   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 421   fh1PtHard->Fill(ptHard,eventW);
 422   fh1PtHard_NoW->Fill(ptHard,1);
 423   fh1PtHard_Trials->Fill(ptHard,nTrials);
 424
 425   // If we set a second branch for the input jets fetch this
 426   if(fBranchGen.Length()>0){
 427     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 428     if(aodGenJets){
 429       Int_t iCount = 0;
 430       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 431         if(iCount>=kMaxJets)continue;
 432         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 433         if(!tmp)continue;
 434         if(fLimitGenJetEta){
 435           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 436              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 437         }
 438         genJets[iCount] = *tmp;
 439         iCount++;
 440       }
 441       nGenJets = iCount;
 442     }
 443     else{
 444       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 445     }
 446   }
 447
 448   fh1NGenJets->Fill(nGenJets);
 449   // We do not want to exceed the maximum jet number
 450   nGenJets = TMath::Min(nGenJets,kMaxJets);
 451
 452   // Fetch the reconstructed jets...
 453
 454
 455   nRecJets = aodRecJets->GetEntries();
 456   fh1NRecJets->Fill(nRecJets);
 457   nRecJets = TMath::Min(nRecJets,kMaxJets);
 458
 459   for(int ir = 0;ir < nRecJets;++ir){
 460     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 461     if(!tmp)continue;
 462     recJets[ir] = *tmp;
 463   }
 464
 465
 466   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 467   // Relate the jets
 468   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 469   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 470
 471   for(int i = 0;i<kMaxJets;++i){
 472     iGenIndex[i] = iRecIndex[i] = -1;
 473   }
 474
 475
 476   GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 477                  iGenIndex,iRecIndex,fDebug);
 478   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 479
 480   if(fDebug){
 481     for(int i = 0;i<kMaxJets;++i){
 482       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 483       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 484     }
 485   }
 486
 487   // loop over reconstructed jets
 488   for(int ir = 0;ir < nRecJets;++ir){
 489     Double_t ptRec = recJets[ir].Pt();
 490     Double_t phiRec = recJets[ir].Phi();
 491     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 492     Double_t etaRec = recJets[ir].Eta();
 493
 494     fh1E[ir]->Fill(recJets[ir].E(),eventW);
 495     fh1PtRecIn[ir]->Fill(ptRec,eventW);
 496     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 497     // Fill Correlation
 498     Int_t ig = iGenIndex[ir];
 499     if(ig>=0&&ig<nGenJets){
 500       fh1PtRecOut[ir]->Fill(ptRec,eventW);
 501       Double_t ptGen = genJets[ig].Pt();
 502       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
 503       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
 504       fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
 505     }
 506     else{
 507       fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 508     }
 509   }// loop over reconstructed jets
 510
 511
 512   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 513   for(int ig = 0;ig < nGenJets;++ig){
 514     Double_t ptGen = genJets[ig].Pt();
 515     // Fill Correlation
 516     Double_t phiGen = genJets[ig].Phi();
 517     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 518     Double_t etaGen = genJets[ig].Eta();
 519     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 520     fh1PtGenIn[ig]->Fill(ptGen,eventW);
 521     Int_t ir = iRecIndex[ig];
 522     if(ir>=0&&ir<nRecJets){
 523       fh1PtGenOut[ig]->Fill(ptGen,eventW);
 524     }
 525     else{
 526       fh3GenEtaPhiPt_NoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 527     }
 528   }// loop over reconstructed jets
 529
 530   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 531   PostData(1, fHistList);
 532 }
 533
 534 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
 535 {
 536 // Terminate analysis
 537 //
 538     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
 539 }
 540
 541
 542 void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
 543                                                 AliAODJet *recJets,Int_t &nRecJets,
 544                                                 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
 545
 546   //
 547   // Relate the two input jet Arrays
 548   //
 549
 550   //
 551   // The association has to be unique
 552   // So check in two directions
 553   // find the closest rec to a gen
 554   // and check if there is no other rec which is closer
 555   // Caveat: Close low energy/split jets may disturb this correlation
 556
 557   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
 558   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
 559   // in the end we have something like this
 560   //    1   2   3
 561   // ------------
 562   // a| 3   2   0
 563   // b| 0   1   0
 564   // c| 0   0   3
 565   // d| 0   0   1
 566   // e| 0   0   1
 567   // Topology
 568   //   1     2
 569   //     a         b
 570   //
 571   //  d      c
 572   //        3     e
 573   // Only entries with "3" match from both sides
 574
 575   const int kMode = 3;
 576
 577   Int_t iFlag[kMaxJets][kMaxJets];
 578
 579
 580
 581   for(int i = 0;i < kMaxJets;++i){
 582     iRecIndex[i] = -1;
 583     iGenIndex[i] = -1;
 584     for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
 585   }
 586
 587   if(nRecJets==0)return;
 588   if(nGenJets==0)return;
 589
 590   const Float_t maxDist = 1.0;
 591   // find the closest distance to the generated
 592   for(int ig = 0;ig<nGenJets;++ig){
 593     Float_t dist = maxDist;
 594     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());
 595     for(int ir = 0;ir<nRecJets;++ir){
 596       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
 597       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());
 598       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
 599       if(dR<dist){
 600         iRecIndex[ig] = ir;
 601         dist = dR;
 602       }
 603     }
 604     if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
 605     // reset...
 606     iRecIndex[ig] = -1;
 607   }
 608   // other way around
 609   for(int ir = 0;ir<nRecJets;++ir){
 610     Float_t dist = maxDist;
 611     for(int ig = 0;ig<nGenJets;++ig){
 612       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
 613       if(dR<dist){
 614         iGenIndex[ir] = ig;
 615         dist = dR;
 616       }
 617     }
 618     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
 619     // reset...
 620     iGenIndex[ir] = -1;
 621   }
 622
 623   // check for "true" correlations
 624
 625   if(iDebug>1)Printf(">>>>>> Matrix");
 626
 627   for(int ig = 0;ig<nGenJets;++ig){
 628     for(int ir = 0;ir<nRecJets;++ir){
 629       // Print
 630       if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
 631
 632       if(kMode==3){
 633         // we have a uniqie correlation
 634         if(iFlag[ig][ir]==3){
 635           iGenIndex[ir] = ig;
 636           iRecIndex[ig] = ir;
 637         }
 638       }
 639       else{
 640         // we just take the correlation from on side
 641         if((iFlag[ig][ir]&2)==2){
 642           iGenIndex[ir] = ig;
 643         }
 644         if((iFlag[ig][ir]&1)==1){
 645           iRecIndex[ig] = ir;
 646         }
 647       }
 648     }
 649     if(iDebug>1)printf("\n");
 650   }
 651 }
 652
 653