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