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