PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.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 <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TList.h>
  35 #include <TLorentzVector.h>
  36 #include <TClonesArray.h>
  37 #include  "TDatabasePDG.h"
  38
  39 #include "AliAnalysisTaskJetSpectrum2.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliJetFinder.h"
  42 #include "AliJetHeader.h"
  43 #include "AliJetReader.h"
  44 #include "AliJetReaderHeader.h"
  45 #include "AliUA1JetHeaderV1.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODTrack.h"
  50 #include "AliAODJet.h"
  51 #include "AliAODMCParticle.h"
  52 #include "AliMCEventHandler.h"
  53 #include "AliMCEvent.h"
  54 #include "AliStack.h"
  55 #include "AliGenPythiaEventHeader.h"
  56 #include "AliJetKineReaderHeader.h"
  57 #include "AliGenCocktailEventHeader.h"
  58 #include "AliInputEventHandler.h"
  59
  60
  61 #include "AliAnalysisHelperJetTasks.h"
  62
  63 ClassImp(AliAnalysisTaskJetSpectrum2)
  64
  65 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
  66   fJetHeaderRec(0x0),
  67   fJetHeaderGen(0x0),
  68   fAOD(0x0),
  69   fhnCorrelation(0x0),
  70   fBranchRec("jets"),
  71   fBranchGen(""),
  72   fUseAODInput(kFALSE),
  73   fUseGlobalSelection(kFALSE),
  74   fUseExternalWeightOnly(kFALSE),
  75   fLimitGenJetEta(kFALSE),
  76   fFilterMask(0),
  77   fAnalysisType(0),
  78   fTrackTypeRec(kTrackUndef),
  79   fTrackTypeGen(kTrackUndef),
  80   fAvgTrials(1),
  81   fExternalWeight(1),
  82   fRecEtaWindow(0.5),
  83   fh1Xsec(0x0),
  84   fh1Trials(0x0),
  85   fh1PtHard(0x0),
  86   fh1PtHardNoW(0x0),
  87   fh1PtHardTrials(0x0),
  88   fh1NGenJets(0x0),
  89   fh1NRecJets(0x0),
  90   fHistList(0x0)
  91 {
  92   for(int i = 0;i < kMaxStep*2;++i){
  93     fhnJetContainer[i] = 0;
  94   }
  95   for(int i = 0;i < kMaxJets;++i){
  96     fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
  97     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
  98   }
  99
 100 }
 101
 102 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 103   AliAnalysisTaskSE(name),
 104   fJetHeaderRec(0x0),
 105   fJetHeaderGen(0x0),
 106   fAOD(0x0),
 107   fhnCorrelation(0x0),
 108   fBranchRec("jets"),
 109   fBranchGen(""),
 110   fUseAODInput(kFALSE),
 111   fUseGlobalSelection(kFALSE),
 112   fUseExternalWeightOnly(kFALSE),
 113   fLimitGenJetEta(kFALSE),
 114   fFilterMask(0),
 115   fAnalysisType(0),
 116   fTrackTypeRec(kTrackUndef),
 117   fTrackTypeGen(kTrackUndef),
 118   fAvgTrials(1),
 119   fExternalWeight(1),
 120   fRecEtaWindow(0.5),
 121   fh1Xsec(0x0),
 122   fh1Trials(0x0),
 123   fh1PtHard(0x0),
 124   fh1PtHardNoW(0x0),
 125   fh1PtHardTrials(0x0),
 126   fh1NGenJets(0x0),
 127   fh1NRecJets(0x0),
 128   fHistList(0x0)
 129 {
 130
 131   for(int i = 0;i < kMaxStep*2;++i){
 132     fhnJetContainer[i] = 0;
 133   }
 134   for(int i = 0;i < kMaxJets;++i){
 135     fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
 136     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 137   }
 138   DefineOutput(1, TList::Class());
 139 }
 140
 141
 142
 143 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 144 {
 145   //
 146   // Implemented Notify() to read the cross sections
 147   // and number of trials from pyxsec.root
 148   //
 149
 150   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 151   Float_t xsection = 0;
 152   Float_t ftrials  = 1;
 153
 154   fAvgTrials = 1;
 155   if(tree){
 156     TFile *curfile = tree->GetCurrentFile();
 157     if (!curfile) {
 158       Error("Notify","No current file");
 159       return kFALSE;
 160     }
 161     if(!fh1Xsec||!fh1Trials){
 162       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 163       return kFALSE;
 164     }
 165     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 166     fh1Xsec->Fill("<#sigma>",xsection);
 167     // construct a poor man average trials
 168     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 169     if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries;
 170   }
 171   return kTRUE;
 172 }
 173
 174 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 175 {
 176
 177   //
 178   // Create the output container
 179   //
 180
 181
 182   // Connect the AOD
 183
 184
 185   MakeJetContainer();
 186
 187
 188   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 189
 190   OpenFile(1);
 191   if(!fHistList)fHistList = new TList();
 192
 193   Bool_t oldStatus = TH1::AddDirectoryStatus();
 194   TH1::AddDirectory(kFALSE);
 195
 196   //
 197   //  Histogram
 198
 199   const Int_t nBinPt = 100;
 200   Double_t binLimitsPt[nBinPt+1];
 201   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 202     if(iPt == 0){
 203       binLimitsPt[iPt] = 0.0;
 204     }
 205     else {// 1.0
 206       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 207     }
 208   }
 209
 210   const Int_t nBinPhi = 30;
 211   Double_t binLimitsPhi[nBinPhi+1];
 212   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 213     if(iPhi==0){
 214       binLimitsPhi[iPhi] = 0;
 215     }
 216     else{
 217       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 218     }
 219   }
 220
 221   const Int_t nBinFrag = 25;
 222
 223
 224   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 225   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 226
 227   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 228   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 229
 230   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 231   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 232   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 233
 234   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 235   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 236
 237   //
 238   for(int ij = 0;ij < kMaxJets;++ij){
 239     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 240     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 241
 242     fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
 243                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 244     fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 245                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 246
 247     fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
 248                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 249     fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 250                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 251   }
 252
 253
 254   const Int_t saveLevel = 3; // large save level more histos
 255   if(saveLevel>0){
 256     fHistList->Add(fh1Xsec);
 257     fHistList->Add(fh1Trials);
 258     fHistList->Add(fh1PtHard);
 259     fHistList->Add(fh1PtHardNoW);
 260     fHistList->Add(fh1PtHardTrials);
 261     fHistList->Add(fh1NGenJets);
 262     fHistList->Add(fh1NRecJets);
 263     for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 264     for(int ij = 0;ij<kMaxJets;++ij){
 265       fHistList->Add( fh1PtRecIn[ij]);
 266       fHistList->Add( fh1PtGenIn[ij]);
 267       fHistList->Add( fh2FragRec[ij]);
 268       fHistList->Add( fh2FragLnRec[ij]);
 269       fHistList->Add( fh2FragGen[ij]);
 270       fHistList->Add( fh2FragLnGen[ij]);
 271     }
 272     fHistList->Add(fhnCorrelation);
 273   }
 274
 275   // =========== Switch on Sumw2 for all histos ===========
 276   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 277     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 278     if (h1){
 279       h1->Sumw2();
 280       continue;
 281     }
 282     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 283     if(hn)hn->Sumw2();
 284   }
 285   TH1::AddDirectory(oldStatus);
 286 }
 287
 288 void AliAnalysisTaskJetSpectrum2::Init()
 289 {
 290   //
 291   // Initialization
 292   //
 293
 294   Printf(">>> AnalysisTaskJetSpectrum2::Init() debug level %d\n",fDebug);
 295   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 296
 297 }
 298
 299 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 300 {
 301
 302   if(! AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
 303     // no selection by the service task, we continue
 304     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 305     PostData(1, fHistList);
 306     return;
 307   }
 308
 309   //
 310   // Execute analysis for current event
 311   //
 312   AliESDEvent *fESD = 0;
 313   if(fUseAODInput){
 314     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 315     if(!fAOD){
 316       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 317       return;
 318     }
 319     // fethc the header
 320   }
 321   else{
 322     //  assume that the AOD is in the general output...
 323     fAOD  = AODEvent();
 324     if(!fAOD){
 325       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 326       return;
 327     }
 328     if(fDebug>0){
 329       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 330     }
 331
 332
 333   }
 334
 335
 336
 337
 338   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 339
 340
 341   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 342
 343   if(!aodH){
 344     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 345     return;
 346   }
 347
 348   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 349   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 350   if(!aodRecJets){
 351     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 352     return;
 353   }
 354
 355   // ==== General variables needed
 356
 357
 358   // We use statice array, not to fragment the memory
 359   AliAODJet genJets[kMaxJets];
 360   Int_t nGenJets = 0;
 361   AliAODJet recJets[kMaxJets];
 362   Int_t nRecJets = 0;
 363   ///////////////////////////
 364
 365
 366   Double_t eventW = 1;
 367   Double_t ptHard = 0;
 368   Double_t nTrials = 1; // Trials for MC trigger
 369
 370   if(fUseExternalWeightOnly){
 371     eventW = fExternalWeight;
 372   }
 373
 374   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 375
 376   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 377   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 378     // this is the part we only use when we have MC information
 379     AliMCEvent* mcEvent = MCEvent();
 380     //    AliStack *pStack = 0;
 381     if(!mcEvent){
 382       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 383       return;
 384     }
 385     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 386     Int_t iCount = 0;
 387     if(pythiaGenHeader){
 388       nTrials = pythiaGenHeader->Trials();
 389       ptHard  = pythiaGenHeader->GetPtHard();
 390       int iProcessType = pythiaGenHeader->ProcessType();
 391       // 11 f+f -> f+f
 392       // 12 f+barf -> f+barf
 393       // 13 f+barf -> g+g
 394       // 28 f+g -> f+g
 395       // 53 g+g -> f+barf
 396       // 68 g+g -> g+g
 397       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 398       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 399
 400       // fetch the pythia generated jets only to be used here
 401       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 402       AliAODJet pythiaGenJets[kMaxJets];
 403       for(int ip = 0;ip < nPythiaGenJets;++ip){
 404         if(iCount>=kMaxJets)continue;
 405         Float_t p[4];
 406         pythiaGenHeader->TriggerJet(ip,p);
 407         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 408
 409         if(fLimitGenJetEta){
 410           if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 411              pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 412         }
 413
 414
 415         if(fBranchGen.Length()==0){
 416           // if we have MC particles and we do not read from the aod branch
 417           // use the pythia jets
 418           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 419         }
 420         iCount++;
 421       }
 422     }
 423     if(fBranchGen.Length()==0)nGenJets = iCount;
 424   }// (fAnalysisType&kMCESD)==kMCESD)
 425
 426
 427   // we simply fetch the tracks/mc particles as a list of AliVParticles
 428
 429   TList recParticles;
 430   TList genParticles;
 431
 432   GetListOfTracks(&recParticles,fTrackTypeRec);
 433   GetListOfTracks(&genParticles,fTrackTypeGen);
 434
 435   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 436   fh1PtHard->Fill(ptHard,eventW);
 437   fh1PtHardNoW->Fill(ptHard,1);
 438   fh1PtHardTrials->Fill(ptHard,nTrials);
 439
 440   // If we set a second branch for the input jets fetch this
 441   if(fBranchGen.Length()>0){
 442     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 443     if(aodGenJets){
 444       Int_t iCount = 0;
 445       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 446         if(iCount>=kMaxJets)continue;
 447         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 448         if(!tmp)continue;
 449         if(fLimitGenJetEta){
 450           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 451              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 452         }
 453         genJets[iCount] = *tmp;
 454         iCount++;
 455       }
 456       nGenJets = iCount;
 457     }
 458     else{
 459       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 460       fAOD->Print();
 461     }
 462   }
 463
 464   fh1NGenJets->Fill(nGenJets);
 465   // We do not want to exceed the maximum jet number
 466   nGenJets = TMath::Min(nGenJets,kMaxJets);
 467
 468   // Fetch the reconstructed jets...
 469
 470   nRecJets = aodRecJets->GetEntries();
 471
 472
 473
 474
 475   fh1NRecJets->Fill(nRecJets);
 476   nRecJets = TMath::Min(nRecJets,kMaxJets);
 477
 478   for(int ir = 0;ir < nRecJets;++ir){
 479     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 480     if(!tmp)continue;
 481     recJets[ir] = *tmp;
 482   }
 483
 484
 485   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 486   // Relate the jets
 487   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 488   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 489
 490   for(int i = 0;i<kMaxJets;++i){
 491     iGenIndex[i] = iRecIndex[i] = -1;
 492   }
 493
 494
 495   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 496                                             iGenIndex,iRecIndex,fDebug);
 497   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 498
 499   if(fDebug){
 500     for(int i = 0;i<kMaxJets;++i){
 501       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 502       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 503     }
 504   }
 505
 506
 507
 508
 509   Double_t container[6];
 510
 511   // loop over generated jets
 512
 513   // radius; tmp, get from the jet header itself
 514   // at some point, how todeal woht FastJet on MC?
 515   Float_t radiusGen =0.4;
 516   Float_t radiusRec =0.4;
 517
 518   for(int ig = 0;ig < nGenJets;++ig){
 519     Double_t ptGen = genJets[ig].Pt();
 520     Double_t phiGen = genJets[ig].Phi();
 521     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 522     Double_t etaGen = genJets[ig].Eta();
 523
 524     container[3] = ptGen;
 525     container[4] = etaGen;
 526     container[5] = phiGen;
 527     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
 528     Int_t ir = iRecIndex[ig];
 529
 530     if(TMath::Abs(etaGen)<fRecEtaWindow){
 531       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
 532       fh1PtGenIn[ig]->Fill(ptGen,eventW);
 533       // fill the fragmentation function
 534       for(int it = 0;it<genParticles.GetEntries();++it){
 535         AliVParticle *part = (AliVParticle*)genParticles.At(it);
 536         if(genJets[ig].DeltaR(part)<radiusGen){
 537           Float_t z = part->Pt()/ptGen;
 538           Float_t lnz =  -1.*TMath::Log(z);
 539           fh2FragGen[ig]->Fill(z,ptGen,eventW);
 540           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
 541         }
 542       }
 543       if(ir>=0&&ir<nRecJets){
 544         fhnJetContainer[kStep3]->Fill(&container[3],eventW);
 545       }
 546     }
 547
 548     if(ir>=0&&ir<nRecJets){
 549       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
 550       Double_t etaRec = recJets[ir].Eta();
 551       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
 552     }
 553   }// loop over generated jets
 554
 555
 556
 557
 558   // loop over reconstructed jets
 559   for(int ir = 0;ir < nRecJets;++ir){
 560     Double_t etaRec = recJets[ir].Eta();
 561     Double_t ptRec = recJets[ir].Pt();
 562     Double_t phiRec = recJets[ir].Phi();
 563     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 564     container[0] = ptRec;
 565     container[1] = etaRec;
 566     container[2] = phiRec;
 567
 568     if(ptRec>20.&&fDebug>0){
 569       // need to cast to int, otherwise the printf overwrites
 570       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
 571       fAOD->GetHeader()->Print();
 572       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
 573         AliAODTrack *tr = fAOD->GetTrack(it);
 574         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 575         tr->Print();
 576         tr->Dump();
 577         if(fESD){
 578           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
 579           esdTr->Print("");
 580           esdTr->Dump();
 581         }
 582       }
 583     }
 584
 585
 586     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
 587     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 588     if(TMath::Abs(etaRec)<fRecEtaWindow){
 589       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
 590       fh1PtRecIn[ir]->Fill(ptRec,eventW);
 591       // fill the fragmentation function
 592       for(int it = 0;it<recParticles.GetEntries();++it){
 593         AliVParticle *part = (AliVParticle*)recParticles.At(it);
 594         if(recJets[ir].DeltaR(part)<radiusRec){
 595           Float_t z = part->Pt()/ptRec;
 596           Float_t lnz =  -1.*TMath::Log(z);
 597           fh2FragRec[ir]->Fill(z,ptRec,eventW);
 598           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
 599         }
 600       }
 601     }
 602
 603
 604     // Fill Correlation
 605     Int_t ig = iGenIndex[ir];
 606     if(ig>=0 && ig<nGenJets){
 607       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
 608       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 609       Double_t ptGen  = genJets[ig].Pt();
 610       Double_t phiGen = genJets[ig].Phi();
 611       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 612       Double_t etaGen = genJets[ig].Eta();
 613
 614       container[3] = ptGen;
 615       container[4] = etaGen;
 616       container[5] = phiGen;
 617
 618       //
 619       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
 620       //
 621
 622       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
 623       if(TMath::Abs(etaRec)<fRecEtaWindow){
 624         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
 625         fhnCorrelation->Fill(container);
 626       }// if etarec in window
 627
 628     }
 629     ////////////////////////////////////////////////////
 630     else{
 631
 632     }
 633   }// loop over reconstructed jets
 634
 635
 636   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 637   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 638   PostData(1, fHistList);
 639 }
 640
 641
 642 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
 643   //
 644   // Create the particle container for the correction framework manager and
 645   // link it
 646   //
 647   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 648   const Double_t kPtmin = 10.0, kPtmax = 260.; // we do not want to have empty bins at the beginning...
 649   const Double_t kEtamin = -3.0, kEtamax = 3.0;
 650   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 651
 652   // can we neglect migration in eta and phi?
 653   // phi should be no problem since we cover full phi and are phi symmetric
 654   // eta migration is more difficult  due to needed acceptance correction
 655   // in limited eta range
 656
 657
 658   //arrays for the number of bins in each dimension
 659   Int_t iBin[kNvar];
 660   iBin[0] = 100; //bins in pt
 661   iBin[1] =  1; //bins in eta
 662   iBin[2] = 1; // bins in phi
 663
 664   //arrays for lower bounds :
 665   Double_t* binEdges[kNvar];
 666   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 667     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 668
 669   //values for bin lower bounds
 670   //  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
 671   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/iBin[0]*(Double_t)i;
 672   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 673   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 674
 675
 676   for(int i = 0;i < kMaxStep*2;++i){
 677     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
 678     for (int k=0; k<kNvar; k++) {
 679       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
 680     }
 681   }
 682   //create correlation matrix for unfolding
 683   Int_t thnDim[2*kNvar];
 684   for (int k=0; k<kNvar; k++) {
 685     //first half  : reconstructed
 686     //second half : MC
 687     thnDim[k]      = iBin[k];
 688     thnDim[k+kNvar] = iBin[k];
 689   }
 690
 691   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
 692   for (int k=0; k<kNvar; k++) {
 693     fhnCorrelation->SetBinEdges(k,binEdges[k]);
 694     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
 695   }
 696   fhnCorrelation->Sumw2();
 697
 698   // Add a histogram for Fake jets
 699   //  thnDim[3] = AliPID::kSPECIES;
 700   //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
 701   // for(Int_t idim = 0; idim < kNvar; idim++)
 702   //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
 703 }
 704
 705 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
 706 {
 707 // Terminate analysis
 708 //
 709     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
 710 }
 711
 712
 713 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
 714
 715
 716   Int_t iCount = 0;
 717   if(type==kTrackAODIn||type==kTrackAODOut){
 718     AliAODEvent *aod = 0;
 719     if(type==kTrackAODIn)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 720     else if (type==kTrackAODOut) aod = AODEvent();
 721     if(!aod){
 722       return iCount;
 723     }
 724     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 725       AliAODTrack *tr = aod->GetTrack(it);
 726       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 727       list->Add(tr);
 728       iCount++;
 729     }
 730   }
 731   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 732     AliMCEvent* mcEvent = MCEvent();
 733     if(!mcEvent)return iCount;
 734     // we want to have alivpartilces so use get track
 735     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 736       if(!mcEvent->IsPhysicalPrimary(it))continue;
 737       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 738       if(type == kTrackKineAll){
 739         list->Add(part);
 740         iCount++;
 741       }
 742       else if(type == kTrackKineCharged){
 743         if(part->Particle()->GetPDG()->Charge()==0)continue;
 744         list->Add(part);
 745         iCount++;
 746       }
 747     }
 748   }
 749   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll ) {
 750     AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
 751     if(!aod) aod = AODEvent();
 752     if(!aod)return iCount;
 753     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 754     if(!tca)return iCount;
 755     for(int it = 0;it < tca->GetEntriesFast();++it){
 756       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 757       if(!part->IsPhysicalPrimary())continue;
 758       if(type == kTrackAODMCAll){
 759         list->Add(part);
 760         iCount++;
 761       }
 762       else if (type == kTrackAODMCCharged){
 763         if(part->Charge()==0)continue;
 764         list->Add(part);
 765         iCount++;
 766       }
 767     }
 768   }// AODMCparticle
 769   return iCount;
 770
 771 }