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