PWGJE/AliAnalysisTaskJetSpectrum2.cxx

   1 // **************************************
   2 // 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 <TRandom3.h>
  26 #include <TSystem.h>
  27 #include <TInterpreter.h>
  28 #include <TChain.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TProfile2D.h>
  36 #include <TList.h>
  37 #include <TLorentzVector.h>
  38 #include <TClonesArray.h>
  39 #include <TRefArray.h>
  40 #include  "TDatabasePDG.h"
  41
  42 #include "AliAnalysisTaskJetSpectrum2.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliJetFinder.h"
  45 #include "AliTHn.h"
  46 #include "AliJetHeader.h"
  47 #include "AliJetReader.h"
  48 #include "AliJetReaderHeader.h"
  49 #include "AliUA1JetHeaderV1.h"
  50 #include "AliESDEvent.h"
  51 #include "AliAODEvent.h"
  52 #include "AliAODHandler.h"
  53 #include "AliAODTrack.h"
  54 #include "AliAODJet.h"
  55 #include "AliAODJetEventBackground.h"
  56 #include "AliAODMCParticle.h"
  57 #include "AliMCEventHandler.h"
  58 #include "AliMCEvent.h"
  59 #include "AliStack.h"
  60 #include "AliGenPythiaEventHeader.h"
  61 #include "AliJetKineReaderHeader.h"
  62 #include "AliGenCocktailEventHeader.h"
  63 #include "AliInputEventHandler.h"
  64 #include "AliCFContainer.h"
  65
  66 #include "AliAnalysisHelperJetTasks.h"
  67
  68 ClassImp(AliAnalysisTaskJetSpectrum2)
  69
  70 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
  71 AliAnalysisTaskSE(),
  72   fJetHeaderRec(0x0),
  73   fJetHeaderGen(0x0),
  74   fAODIn(0x0),
  75   fAODOut(0x0),
  76   fAODExtension(0x0),
  77   fhnJetContainer(0x0),
  78   fhnCorrelation(0x0),
  79   fhnEvent(0x0),
  80   f1PtScale(0x0),
  81   fBranchRec("jets"),
  82   fBranchGen(""),
  83   fBranchBkgRec(""),
  84   fBranchBkgGen(""),
  85   fNonStdFile(""),
  86   fRandomizer(0x0),
  87   fUseAODJetInput(kFALSE),
  88   fUseAODTrackInput(kFALSE),
  89   fUseAODMCInput(kFALSE),
  90   fUseGlobalSelection(kFALSE),
  91   fUseExternalWeightOnly(kFALSE),
  92   fLimitGenJetEta(kFALSE),
  93   fDoMatching(kFALSE),
  94   fNMatchJets(5),
  95   fNRPBins(3),
  96   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
  97   fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
  98   fFilterMask(0),
  99   fEventSelectionMask(0),
 100   fNTrigger(0),
 101   fTriggerBit(0x0),
 102   fNAcceptance(0),
 103   fNBinsLeadingTrackPt(10),
 104   fNBinsMult(20),
 105   fAnalysisType(0),
 106   fTrackTypeRec(kTrackUndef),
 107   fTrackTypeGen(kTrackUndef),
 108   fEventClass(0),
 109   fRPMethod(0),
 110   fAvgTrials(1),
 111   fExternalWeight(1),
 112   fJetRecEtaWindow(0.5),
 113   fTrackRecEtaWindow(0.5),
 114   fMinJetPt(0),
 115   fMinTrackPt(0.15),
 116   fDeltaPhiWindow(90./180.*TMath::Pi()),
 117   fAcceptancePhiMin(0x0),
 118   fAcceptancePhiMax(0x0),
 119   fAcceptanceEtaMin(0x0),
 120   fAcceptanceEtaMax(0x0),
 121   fCentrality(100),
 122   fRPAngle(0),
 123   fMultRec(0),
 124   fMultGen(0),
 125   fTriggerName(0x0),
 126   fh1Xsec(0x0),
 127   fh1Trials(0x0),
 128   fh1AvgTrials(0x0),
 129   fh1PtHard(0x0),
 130   fh1PtHardNoW(0x0),
 131   fh1PtHardTrials(0x0),
 132   fh1ZVtx(0x0),
 133   fh1RP(0x0),
 134   fh1Centrality(0x0),
 135   fh1TmpRho(0x0),
 136   fh2MultRec(0x0),
 137   fh2MultGen(0x0),
 138   fh2RPCentrality(0x0),
 139   fh2PtFGen(0x0),
 140   fh2RelPtFGen(0x0),
 141   fh3RelPtFGenLeadTrkPt(0x0),
 142   fHistList(0x0)
 143 {
 144
 145   for(int ij = 0;ij <kJetTypes;++ij){
 146     fFlagJetType[ij] = 1; // default = on
 147     fh1NJets[ij] = 0;
 148     fh1SumPtTrack[ij] = 0;
 149     fh1PtJetsIn[ij] = 0;
 150     fh1PtJetsInRej[ij] = 0;
 151     fh1PtJetsInBest[ij] = 0;
 152     fh1PtTracksIn[ij] = 0;
 153     fh1PtTracksInLow[ij] = 0;
 154     fh2NJetsPt[ij]  = 0;
 155     fh2NTracksPt[ij]  = 0;
 156     fp2MultRPPhiTrackPt[ij] = 0;
 157     fp2CentRPPhiTrackPt[ij] = 0;
 158     fhnJetPt[ij] = 0;
 159     fhnJetPtBest[ij] = 0;
 160     fhnJetPtRej[ij] = 0;
 161     fhnJetPtQA[ij] = 0;
 162     fhnTrackPt[ij] = 0;
 163     fhnTrackPtQA[ij] = 0;
 164     for(int i = 0;i <= kMaxJets;++i){
 165       fh2LTrackPtJetPt[ij][i] = 0;
 166       fh1PtIn[ij][i] = 0;
 167     }
 168
 169     fh1DijetMinv[ij] = 0;
 170     fh2DijetDeltaPhiPt[ij] = 0;
 171     fh2DijetAsymPt[ij] = 0;
 172     fh2DijetPt2vsPt1[ij] = 0;
 173     fh2DijetDifvsSum[ij] = 0;
 174
 175   }
 176 }
 177
 178 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 179   AliAnalysisTaskSE(name),
 180   fJetHeaderRec(0x0),
 181   fJetHeaderGen(0x0),
 182   fAODIn(0x0),
 183   fAODOut(0x0),
 184   fAODExtension(0x0),
 185   fhnJetContainer(0x0),
 186   fhnCorrelation(0x0),
 187   fhnEvent(0x0),
 188   f1PtScale(0x0),
 189   fBranchRec("jets"),
 190   fBranchGen(""),
 191   fBranchBkgRec(""),
 192   fBranchBkgGen(""),
 193   fNonStdFile(""),
 194   fRandomizer(0x0),
 195   fUseAODJetInput(kFALSE),
 196   fUseAODTrackInput(kFALSE),
 197   fUseAODMCInput(kFALSE),
 198   fUseGlobalSelection(kFALSE),
 199   fUseExternalWeightOnly(kFALSE),
 200   fLimitGenJetEta(kFALSE),
 201   fDoMatching(kFALSE),
 202   fNMatchJets(5),
 203   fNRPBins(3),
 204   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 205   fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
 206   fFilterMask(0),
 207   fEventSelectionMask(0),
 208   fNTrigger(0),
 209   fTriggerBit(0x0),
 210   fNAcceptance(0),
 211   fNBinsLeadingTrackPt(10),
 212   fNBinsMult(20),
 213   fAnalysisType(0),
 214   fTrackTypeRec(kTrackUndef),
 215   fTrackTypeGen(kTrackUndef),
 216   fEventClass(0),
 217   fRPMethod(0),
 218   fAvgTrials(1),
 219   fExternalWeight(1),
 220   fJetRecEtaWindow(0.5),
 221   fTrackRecEtaWindow(0.5),
 222   fMinJetPt(0),
 223   fMinTrackPt(0.15),
 224   fDeltaPhiWindow(90./180.*TMath::Pi()),
 225   fAcceptancePhiMin(0x0),
 226   fAcceptancePhiMax(0x0),
 227   fAcceptanceEtaMin(0x0),
 228   fAcceptanceEtaMax(0x0),
 229   fCentrality(100),
 230   fRPAngle(0),
 231   fMultRec(0),
 232   fMultGen(0),
 233   fTriggerName(0x0),
 234   fh1Xsec(0x0),
 235   fh1Trials(0x0),
 236   fh1AvgTrials(0x0),
 237   fh1PtHard(0x0),
 238   fh1PtHardNoW(0x0),
 239   fh1PtHardTrials(0x0),
 240   fh1ZVtx(0x0),
 241   fh1RP(0x0),
 242   fh1Centrality(0x0),
 243   fh1TmpRho(0x0),
 244   fh2MultRec(0x0),
 245   fh2MultGen(0x0),
 246   fh2RPCentrality(0x0),
 247   fh2PtFGen(0x0),
 248   fh2RelPtFGen(0x0),
 249   fh3RelPtFGenLeadTrkPt(0x0),
 250   fHistList(0x0)
 251 {
 252
 253   for(int ij = 0;ij <kJetTypes;++ij){
 254     fFlagJetType[ij] = 1; // default = on
 255     fh1NJets[ij] = 0;
 256     fh1SumPtTrack[ij] = 0;
 257     fh1PtJetsIn[ij] = 0;
 258     fh1PtJetsInRej[ij] = 0;
 259     fh1PtJetsInBest[ij] = 0;
 260     fh1PtTracksIn[ij] = 0;
 261     fh1PtTracksInLow[ij] = 0;
 262     fh2NJetsPt[ij]  = 0;
 263     fh2NTracksPt[ij]  = 0;
 264     fp2MultRPPhiTrackPt[ij] = 0;
 265     fp2CentRPPhiTrackPt[ij] = 0;
 266     fhnJetPt[ij] = 0;
 267     fhnJetPtBest[ij] = 0;
 268     fhnJetPtRej[ij] = 0;
 269     fhnJetPtQA[ij] = 0;
 270     fhnTrackPt[ij] = 0;
 271     fhnTrackPtQA[ij] = 0;
 272     for(int i = 0;i <= kMaxJets;++i){
 273       fh2LTrackPtJetPt[ij][i] = 0;
 274       fh1PtIn[ij][i] = 0;
 275     }
 276
 277     fh1DijetMinv[ij] = 0;
 278     fh2DijetDeltaPhiPt[ij] = 0;
 279     fh2DijetAsymPt[ij] = 0;
 280     fh2DijetPt2vsPt1[ij] = 0;
 281     fh2DijetDifvsSum[ij] = 0;
 282   }
 283
 284   DefineOutput(1, TList::Class());
 285 }
 286
 287
 288
 289 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 290 {
 291
 292
 293
 294   //
 295   // Implemented Notify() to read the cross sections
 296   // and number of trials from pyxsec.root
 297   //
 298
 299   // Fetch the aod also from the input in,
 300   // have todo it in notify
 301
 302
 303   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 304   //  assume that the AOD is in the general output...
 305   fAODOut  = AODEvent();
 306
 307   if(fNonStdFile.Length()!=0){
 308     // case that we have an AOD extension we need can fetch the jets from the extended output
 309     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 310     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 311     if(!fAODExtension){
 312       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 313     }
 314   }
 315
 316
 317   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 318   Float_t xsection = 0;
 319   Float_t ftrials  = 1;
 320
 321   fAvgTrials = 1;
 322   if(tree){
 323     TFile *curfile = tree->GetCurrentFile();
 324     if (!curfile) {
 325       Error("Notify","No current file");
 326       return kFALSE;
 327     }
 328     if(!fh1Xsec||!fh1Trials){
 329       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 330       return kFALSE;
 331     }
 332     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 333     fh1Xsec->Fill("<#sigma>",xsection);
 334     // construct a poor man average trials
 335     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 336     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 337     fh1Trials->Fill("#sum{ntrials}",ftrials);
 338   }
 339
 340   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 341
 342   return kTRUE;
 343 }
 344
 345 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 346 {
 347
 348
 349   // Connect the AOD
 350
 351   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 352   OpenFile(1);
 353   if(!fHistList)fHistList = new TList();
 354   PostData(1, fHistList); // post data in any case once
 355
 356   if(!fRandomizer)fRandomizer = new TRandom3(0);
 357
 358   fHistList->SetOwner(kTRUE);
 359   Bool_t oldStatus = TH1::AddDirectoryStatus();
 360   TH1::AddDirectory(kFALSE);
 361
 362
 363
 364   // event npsparse cent, mult
 365   const Int_t nBinsSparse0 = 4;
 366   const Int_t nBins0[nBinsSparse0] = {     100, 500,fNTrigger,125};
 367   const Double_t xmin0[nBinsSparse0]  = {    0,   0, -0.5,-2};
 368   const Double_t xmax0[nBinsSparse0]  = {  100,5000,fNTrigger-0.5,248};
 369
 370
 371   fhnEvent = new THnSparseF("fhnEvent",";cent;mult:trigger;#rho",nBinsSparse0,
 372                             nBins0,xmin0,xmax0);
 373   fHistList->Add(fhnEvent);
 374
 375   if(fDoMatching){
 376     MakeJetContainer();
 377     fHistList->Add(fhnCorrelation);
 378     fHistList->Add(fhnJetContainer);
 379   }
 380   //
 381   //  Histogram
 382
 383
 384
 385   const Int_t nBinPt = 150;
 386   Double_t binLimitsPt[nBinPt+1];
 387   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 388     if(iPt == 0){
 389       binLimitsPt[iPt] = -50.0;
 390     }
 391     else {// 1.0
 392       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
 393     }
 394   }
 395   const Int_t nBinPhi = 90;
 396   Double_t binLimitsPhi[nBinPhi+1];
 397   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 398     if(iPhi==0){
 399       binLimitsPhi[iPhi] = 0;
 400     }
 401     else{
 402       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 403     }
 404   }
 405
 406   const Int_t nBinEta = 40;
 407   Double_t binLimitsEta[nBinEta+1];
 408   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 409     if(iEta==0){
 410       binLimitsEta[iEta] = -2.0;
 411     }
 412     else{
 413       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 414     }
 415   }
 416
 417
 418   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 419   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 420   fHistList->Add(fh1Xsec);
 421   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 422   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 423   fHistList->Add(fh1Trials);
 424   fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
 425   fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
 426   fHistList->Add(fh1AvgTrials);
 427   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 428   fHistList->Add(fh1PtHard);
 429   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 430   fHistList->Add(fh1PtHardNoW);
 431   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 432   fHistList->Add(fh1PtHardTrials);
 433
 434   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 435   fHistList->Add(fh1ZVtx);
 436
 437
 438   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
 439   fHistList->Add(fh1RP);
 440
 441   fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
 442   fHistList->Add(fh1Centrality);
 443
 444   fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
 445   fHistList->Add(fh2MultRec);
 446   fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
 447   fHistList->Add(fh2MultGen);
 448
 449
 450   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
 451   fHistList->Add(fh2RPCentrality);
 452
 453   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 454   fHistList->Add(fh2PtFGen);
 455
 456   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
 457   fHistList->Add(fh2RelPtFGen);
 458
 459   const Int_t nBinsLeadingTrackPt = 10;
 460   const Double_t binArrayLeadingTrackPt[nBinsLeadingTrackPt+1] = {0.,1.,2.,3.,4.,5.,6.,8.,10.,12.,200.}; //store pT of leading track in jet
 461
 462   const Int_t nBinDeltaPt = 241;
 463   Double_t binLimitsDeltaPt[nBinDeltaPt+1];
 464   for(Int_t iDeltaPt = 0;iDeltaPt <= nBinDeltaPt;iDeltaPt++){
 465     if(iDeltaPt == 0){
 466       binLimitsDeltaPt[iDeltaPt] = -2.41;
 467     }
 468     else {
 469       binLimitsDeltaPt[iDeltaPt] =  binLimitsDeltaPt[iDeltaPt-1] + 0.02;
 470     }
 471   }
 472
 473   fh3RelPtFGenLeadTrkPt = new TH3F("fh3RelPtFGenLeadTrkPt",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen};p_{T}^{leading track}",nBinPt,binLimitsPt,nBinDeltaPt,binLimitsDeltaPt,nBinsLeadingTrackPt,binArrayLeadingTrackPt);
 474   fHistList->Add(fh3RelPtFGenLeadTrkPt);
 475
 476   for(int ij = 0;ij <kJetTypes;++ij){
 477     TString cAdd = "";
 478     TString cJetBranch = "";
 479     if(ij==kJetRec){
 480       cAdd = "Rec";
 481       cJetBranch = fBranchRec.Data();
 482     }
 483     else if (ij==kJetGen){
 484       cAdd = "Gen";
 485       cJetBranch = fBranchGen.Data();
 486     }
 487     else if (ij==kJetRecFull){
 488       cAdd = "RecFull";
 489       cJetBranch = fBranchRec.Data();
 490     }
 491     else if (ij==kJetGenFull){
 492       cAdd = "GenFull";
 493       cJetBranch = fBranchGen.Data();
 494     }
 495
 496     if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
 497     if(!fFlagJetType[ij])continue;
 498
 499     fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
 500     fHistList->Add(fh1NJets[ij]);
 501
 502     fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 503     fHistList->Add(fh1PtJetsIn[ij]);
 504
 505     fh1PtJetsInRej[ij]  = new TH1F(Form("fh1PtJets%sInRej",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 506     fHistList->Add(fh1PtJetsInRej[ij]);
 507
 508     fh1PtJetsInBest[ij]  = new TH1F(Form("fh1PtJets%sInBest",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 509     fHistList->Add(fh1PtJetsInBest[ij]);
 510
 511     fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 512     fHistList->Add(fh1PtTracksIn[ij]);
 513
 514     fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
 515     fHistList->Add(fh1PtTracksInLow[ij]);
 516
 517     fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),1000,0.,3000.);
 518     fHistList->Add(fh1SumPtTrack[ij]);
 519
 520     fh2NJetsPt[ij]  = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
 521     fHistList->Add(fh2NJetsPt[ij]);
 522
 523     fh2NTracksPt[ij]  = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,5000);
 524     fHistList->Add(fh2NTracksPt[ij]);
 525
 526
 527     fp2MultRPPhiTrackPt[ij] = new TProfile2D(Form("fp2MultRPPhiTrackPt%s",cAdd.Data()),"RP phi vs Number of tracks;# tracks;#Delta#phi_{RP}; <p_{T}>",20,0,4000,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
 528     fHistList->Add(fp2MultRPPhiTrackPt[ij]);
 529     fp2CentRPPhiTrackPt[ij] = new TProfile2D(Form("fp2CentRPPhiTrackPt%s",cAdd.Data()),"RP phi vs cent;# cent;#Delta#phi_{RP}; <p_{T}>",10,0,100,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
 530     fHistList->Add(fp2CentRPPhiTrackPt[ij]);
 531
 532     // Bins:  Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M
 533     const Int_t nBinsSparse1 = 9;
 534     const Int_t nBinsArea = 10;
 535     Int_t nBins1[nBinsSparse1] = {     kMaxJets+1,120, 10,  fNBinsMult,    fNRPBins, nBinsArea,fNTrigger,fNBinsLeadingTrackPt,fNAcceptance+1};
 536     if(cJetBranch.Contains("RandomCone")){
 537       nBins1[1] = 600;
 538       nBins1[5] = 1;
 539     }
 540     const Double_t xmin1[nBinsSparse1]  = {        -0.5,-50,  0,   0,        -0.5, 0.,         -0.5,  0.,           -0.5,};
 541     const Double_t xmax1[nBinsSparse1]  = {kMaxJets+0.5,250,100,4000,fNRPBins-0.5,1.0,fNTrigger-0.5,200.,fNAcceptance+0.5};
 542
 543     const Double_t binArrayArea[nBinsArea+1] = {xmin1[5],0.07,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,xmax1[5]};
 544
 545
 546     fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leading track p_{T};acceptance bin",nBinsSparse1,nBins1,xmin1,xmax1);
 547     fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
 548     if(fNBinsLeadingTrackPt>1) fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
 549     fHistList->Add(fhnJetPt[ij]);
 550
 551
 552     // Bins:  pTJet, cent, trigger,
 553     const Int_t nBinsSparse1b = 3;
 554     Int_t nBins1b[nBinsSparse1b] = {120, 10,fNTrigger};
 555     const Double_t xmin1b[nBinsSparse1b]  = {-50,  0,-0.5};
 556     const Double_t xmax1b[nBinsSparse1b]  = {250,100,fNTrigger-0.5};
 557
 558     fhnJetPtBest[ij] = new THnSparseF(Form("fhnJetPtBest%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
 559     fHistList->Add(fhnJetPtBest[ij]);
 560
 561     fhnJetPtRej[ij] = new THnSparseF(Form("fhnJetPtRej%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
 562     fHistList->Add(fhnJetPtRej[ij]);
 563
 564     // Bins:  Jet number: pTJet, cent, eta, phi, Area, trigger, acceptance, signed pT leading
 565     const Int_t nBinsSparse2 = 9;
 566     Int_t nBins2[nBinsSparse2] = {     kMaxJets+1, 60,   8,  18, 72, 10,fNTrigger,fNAcceptance+1,20};
 567     if(cJetBranch.Contains("RandomCone")){
 568       nBins2[5] = 1;
 569     }
 570     const Double_t xmin2[nBinsSparse2]  = {        -0.5, -50,  0,-0.9,              0,  0.,        -0.5,            -0.5, -100};
 571     const Double_t xmax2[nBinsSparse2]  = {kMaxJets+0.5, 250, 80, 0.9, 2.*TMath::Pi(),1.0,fNTrigger-0.5,fNAcceptance+0.5,  100};
 572     fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading",nBinsSparse2,nBins2,xmin2,xmax2);
 573     fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
 574     fHistList->Add(fhnJetPtQA[ij]);
 575
 576     // Bins:track number  pTtrack, cent, mult, RP.  charge total bins = 224 k
 577     const Int_t nBinsSparse3 = 7;
 578     const Int_t nRPBinsSparse3 = 3;
 579     const Int_t nBins3[nBinsSparse3] = {       2,    100,     10,   1,    nRPBinsSparse3,fNTrigger,2};
 580     const Double_t xmin3[nBinsSparse3]  = { -0.5,      0,   0,      0,        -0.5,-0.5,-1.5};
 581     const Double_t xmax3[nBinsSparse3]  = { 1.5,     200, 100,   4000,nRPBinsSparse3-0.5,fNTrigger-0.5,1.5};
 582
 583     // change the binning of the pT axis:
 584     Double_t *xPt3 = new Double_t[nBins3[1]+1];
 585     xPt3[0] = 0.;
 586     for(int i = 1; i<=nBins3[1];i++){
 587       if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
 588       else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
 589       else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
 590       else if(xPt3[i-1]<15)xPt3[i] = xPt3[i-1] +  1.; // 62 - 67
 591       else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 2.; // 67 - 72
 592       else if(xPt3[i-1]<60)xPt3[i] = xPt3[i-1] + 5; // 72 - 76
 593       else xPt3[i] = xPt3[i-1] + 10; // 76 - 100 = 140
 594     }
 595
 596     fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP;trigger",nBinsSparse3,nBins3,xmin3,xmax3);
 597     fhnTrackPt[ij]->SetBinEdges(1,xPt3);
 598     fHistList->Add(fhnTrackPt[ij]);
 599     delete [] xPt3;
 600
 601     // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M
 602     const Int_t nBinsSparse4 = 6;
 603     const Int_t nBins4[nBinsSparse4] =    {    2, 50,  10,   20, 180,2};
 604     const Double_t xmin4[nBinsSparse4]  = { -0.5,  0,   0, -1.0,  0.,-1.5};
 605     const Double_t xmax4[nBinsSparse4]  = {  1.5,150, 100,  1.0,2.*TMath::Pi(),1.5};
 606
 607     // change the binning ot the pT axis:
 608     Double_t *xPt4 = new Double_t[nBins4[1]+1];
 609     xPt4[0] = 0.;
 610     for(int i = 1; i<=nBins4[1];i++){
 611       if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1;
 612       else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5;
 613       else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] +  1.;
 614       else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] +  2.5;
 615       else xPt4[i] = xPt4[i-1] + 5.;
 616     }
 617     fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi;sign",nBinsSparse4,nBins4,xmin4,xmax4);
 618     fhnTrackPtQA[ij]->SetBinEdges(1,xPt4);
 619     fHistList->Add(fhnTrackPtQA[ij]);
 620     delete [] xPt4;
 621
 622     for(int i = 0;i <= kMaxJets;++i){
 623       fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
 624       fHistList->Add(fh1PtIn[ij][i]);
 625
 626
 627       if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 628       fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
 629                                          Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
 630                                               cAdd.Data()),
 631                                          200,0.,200.,nBinPt,binLimitsPt);
 632       fHistList->Add(fh2LTrackPtJetPt[ij][i]);
 633     }
 634
 635
 636     fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 637     fHistList->Add(fh1DijetMinv[ij]);
 638
 639     fh2DijetDeltaPhiPt[ij]       = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
 640     fHistList->Add(fh2DijetDeltaPhiPt[ij]);
 641
 642     fh2DijetAsymPt[ij]            = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
 643     fHistList->Add(fh2DijetAsymPt[ij]);
 644
 645     fh2DijetPt2vsPt1[ij]          = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
 646     fHistList->Add(fh2DijetPt2vsPt1[ij]);
 647     fh2DijetDifvsSum[ij]         = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
 648     fHistList->Add( fh2DijetDifvsSum[ij]);
 649   }
 650   // =========== Switch on Sumw2 for all histos ===========
 651   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 652     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 653     if (h1){
 654       h1->Sumw2();
 655       continue;
 656     }
 657     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 658     if(hn)hn->Sumw2();
 659   }
 660   TH1::AddDirectory(oldStatus);
 661 }
 662
 663 void AliAnalysisTaskJetSpectrum2::Init()
 664 {
 665   //
 666   // Initialization
 667   //
 668
 669   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 670
 671 }
 672
 673 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
 674
 675
 676   Bool_t selected = kTRUE;
 677
 678   if(fUseGlobalSelection&&fEventSelectionMask==0){
 679     selected = AliAnalysisHelperJetTasks::Selected();
 680   }
 681   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 682     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 683   }
 684
 685   if(fEventClass>0){
 686     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 687   }
 688
 689   if(!selected){
 690     // no selection by the service task, we continue
 691     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 692     PostData(1, fHistList);
 693     return;
 694   }
 695
 696
 697   static AliAODEvent* aod = 0;
 698
 699   // take all other information from the aod we take the tracks from
 700   if(!aod){
 701     if(fUseAODTrackInput)aod = fAODIn;
 702     else aod = fAODOut;
 703   }
 704
 705
 706   //
 707   // Execute analysis for current event
 708   //
 709   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 710   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 711   if(!aodH){
 712     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 713     return;
 714   }
 715
 716   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 717   TClonesArray *aodRecJets = 0;
 718
 719   if(fAODOut&&!aodRecJets){
 720     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
 721   }
 722   if(fAODExtension&&!aodRecJets){
 723     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
 724   }
 725   if(fAODIn&&!aodRecJets){
 726     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
 727   }
 728
 729
 730
 731   if(!aodRecJets){
 732     if(fDebug){
 733
 734       Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 735       if(fAODIn){
 736         Printf("Input AOD >>>>");
 737         fAODIn->Print();
 738       }
 739       if(fAODExtension){
 740         Printf("AOD Extension >>>>");
 741         fAODExtension->Print();
 742       }
 743       if(fAODOut){
 744         Printf("Output AOD >>>>");
 745         fAODOut->Print();
 746       }
 747     }
 748     return;
 749   }
 750
 751   TClonesArray *aodGenJets = 0;
 752   if(fBranchGen.Length()>0){
 753     if(fAODOut&&!aodGenJets){
 754       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
 755     }
 756     if(fAODExtension&&!aodGenJets){
 757       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
 758     }
 759     if(fAODIn&&!aodGenJets){
 760       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
 761     }
 762
 763     if(!aodGenJets){
 764       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
 765       return;
 766     }
 767   }
 768
 769   TClonesArray *aodBackRecJets = 0;
 770   if(fBranchBkgRec.Length()>0){
 771     if(fAODOut&&!aodBackRecJets){
 772       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
 773     }
 774     if(fAODExtension&&!aodBackRecJets){
 775       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
 776     }
 777     if(fAODIn&&!aodBackRecJets){
 778       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
 779     }
 780
 781     if(!aodBackRecJets){
 782       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
 783       return;
 784     }
 785   }
 786
 787
 788   TClonesArray *aodBackGenJets = 0;
 789
 790   if(fBranchBkgGen.Length()>0){
 791     if(fAODOut&&!aodBackGenJets){
 792       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
 793     }
 794     if(fAODExtension&&!aodBackGenJets){
 795       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
 796     }
 797     if(fAODIn&&!aodBackGenJets){
 798       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
 799     }
 800
 801     if(!aodBackGenJets){
 802       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
 803       return;
 804     }
 805   }
 806
 807
 808   // new Scheme
 809   // first fill all the pure  histograms, i.e. full jets
 810   // in case of the correaltion limit to the n-leading jets
 811
 812   // reconstructed
 813
 814
 815   // generated
 816
 817
 818   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
 819
 820
 821
 822   TList genJetsList;         // full acceptance
 823   TList genJetsListCut;      // acceptance cut
 824   TList recJetsList;         // full acceptance
 825   TList recJetsListCut;      // acceptance cut
 826
 827
 828   GetListOfJets(&recJetsList,aodRecJets,0);
 829   GetListOfJets(&recJetsListCut,aodRecJets,1);
 830
 831   if(fBranchGen.Length()>0){
 832     GetListOfJets(&genJetsList,aodGenJets,0);
 833     GetListOfJets(&genJetsListCut,aodGenJets,1);
 834   }
 835
 836   Double_t eventW = 1;
 837   Double_t ptHard = 0;
 838   Double_t nTrials = 1; // Trials for MC trigger
 839   fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
 840
 841   // Getting some global properties
 842   fCentrality = GetCentrality();
 843   if(fCentrality<=0)fCentrality = 0;
 844   fh1Centrality->Fill(fCentrality);
 845
 846
 847   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 848     // this is the part we only use when we have MC information
 849     AliMCEvent* mcEvent = MCEvent();
 850     //    AliStack *pStack = 0;
 851     if(!mcEvent){
 852       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 853       return;
 854     }
 855     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 856     if(pythiaGenHeader){
 857       nTrials = pythiaGenHeader->Trials();
 858       ptHard  = pythiaGenHeader->GetPtHard();
 859       int iProcessType = pythiaGenHeader->ProcessType();
 860       // 11 f+f -> f+f
 861       // 12 f+barf -> f+barf
 862       // 13 f+barf -> g+g
 863       // 28 f+g -> f+g
 864       // 53 g+g -> f+barf
 865       // 68 g+g -> g+g
 866       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 867       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 868     }
 869   }// (fAnalysisType&kMCESD)==kMCESD)
 870
 871
 872   // we simply fetch the tracks/mc particles as a list of AliVParticles
 873
 874   TList recParticles;
 875   TList genParticles;
 876
 877   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 878   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 879   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 880   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 881
 882   //  CalculateReactionPlaneAngle(&recParticles);
 883   fRPAngle = 0;
 884
 885   if(fRPMethod==0)fRPAngle = aod->GetHeader()->GetEventplane();
 886   else if(fRPMethod==1||fRPMethod==2){
 887     fRPAngle = aod->GetHeader()->GetQTheta(fRPMethod);
 888   }
 889   fh1RP->Fill(fRPAngle);
 890   fh2RPCentrality->Fill(fCentrality,fRPAngle);
 891   // Event control and counting ...
 892   // MC
 893   fh1PtHard->Fill(ptHard,eventW);
 894   fh1PtHardNoW->Fill(ptHard,1);
 895   fh1PtHardTrials->Fill(ptHard,nTrials);
 896
 897   // Real
 898   if(aod->GetPrimaryVertex()){// No vtx for pure MC
 899     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
 900   }
 901
 902
 903   Int_t recMult1 = recParticles.GetEntries();
 904   Int_t genMult1 = genParticles.GetEntries();
 905
 906   Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
 907   Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
 908
 909   fh2MultRec->Fill(recMult1,recMult2);
 910   fh2MultGen->Fill(genMult1,genMult2);
 911   fMultRec = recMult1;
 912   if(fMultRec<=0)fMultRec = recMult2;
 913   fMultGen = genMult1;
 914   if(fMultGen<=0)fMultGen = genMult2;
 915
 916   Double_t var0[4] = {0,};
 917   var0[0] = fCentrality;
 918   var0[1] = fMultRec;
 919   for(int it=0;it<fNTrigger;it++){
 920     if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 921       var0[2] = it;
 922       var0[3] = GetRho(recJetsList);
 923       fhnEvent->Fill(var0);
 924     }
 925   }
 926   // the loops for rec and gen should be indentical... pass it to a separate
 927   // function ...
 928   // Jet Loop
 929   // Track Loop
 930   // Jet Jet Loop
 931   // Jet Track loop
 932
 933   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 934   FillJetHistos(recJetsList,recParticles,kJetRecFull);
 935   FillTrackHistos(recParticles,kJetRec);
 936
 937   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 938   FillJetHistos(genJetsList,genParticles,kJetGenFull);
 939   FillTrackHistos(genParticles,kJetGen);
 940
 941   // Here follows the jet matching part
 942   // delegate to a separated method?
 943
 944   if(fDoMatching){
 945     FillMatchHistos(recJetsList,genJetsList);
 946   }
 947
 948   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 949   PostData(1, fHistList);
 950 }
 951
 952 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 953
 954   if(iType>=kJetTypes){
 955     return;
 956   }
 957   if(!fFlagJetType[iType])return;
 958
 959   Int_t refMult = fMultRec;
 960   if(iType==kJetGen||iType==kJetGenFull){
 961     refMult = fMultGen;
 962   }
 963
 964   Int_t nJets = jetsList.GetEntries();
 965   fh1NJets[iType]->Fill(nJets);
 966
 967   if(nJets<=0)return;
 968
 969   Float_t ptOld = 1.E+32;
 970   Float_t pT0 = 0;
 971   Float_t pT1 = 0;
 972   Float_t phi0 = 0;
 973   Float_t phi1 = 0;
 974   Int_t ij0 = -1;
 975   Int_t ij1 = -1;
 976
 977   Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
 978   Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
 979
 980   var1[2] = fCentrality;
 981   var1b[1] = fCentrality;
 982   var1[3] = refMult;
 983
 984
 985
 986   Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
 987   var2[2] = fCentrality;
 988
 989   for(int ij = 0;ij < nJets;ij++){
 990     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
 991     Float_t ptJet = jet->Pt();
 992
 993
 994     if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
 995
 996     var1b[0] = ptJet;
 997     if(jet->Trigger()&fJetTriggerBestMask){
 998       fh1PtJetsInBest[iType]->Fill(ptJet);
 999       for(int it = 0;it <fNTrigger;it++){
1000         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1001           var1b[2] = it;
1002           fhnJetPtBest[iType]->Fill(var1b);
1003         }
1004       }
1005     }
1006     if(jet->Trigger()&fJetTriggerExcludeMask){
1007       fh1PtJetsInRej[iType]->Fill(ptJet);
1008       for(int it = 0;it <fNTrigger;it++){
1009         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1010           var1b[2] = it;
1011           fhnJetPtRej[iType]->Fill(var1b);
1012         }
1013       }
1014       continue;
1015     }
1016     fh1PtJetsIn[iType]->Fill(ptJet);
1017     ptOld = ptJet;
1018
1019     // find the dijets assume sorting and acceptance cut...
1020     if(ij==0){
1021       ij0 = ij;
1022       pT0 = ptJet;
1023       phi0 = jet->Phi();
1024       if(phi0<0)phi0+=TMath::Pi()*2.;
1025     }
1026     else if(ptJet>pT1){
1027       // take only the backward hemisphere??
1028       phi1 = jet->Phi();
1029       if(phi1<0)phi1+=TMath::Pi()*2.;
1030       Float_t dPhi = phi1 - phi0;
1031       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1032       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1033       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1034         ij1 = ij;
1035         pT1 = ptJet;
1036       }
1037     }
1038     // fill jet histos for kmax jets
1039
1040     Float_t phiJet = jet->Phi();
1041     Float_t etaJet = jet->Eta();
1042     if(phiJet<0)phiJet+=TMath::Pi()*2.;
1043     fh1TmpRho->Reset();
1044     if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1045
1046     fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1047     // fill leading jets...
1048     AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1049     //      AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1050     Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1051     Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1052
1053     var1[1] = ptJet;
1054     var1[4] = phiBin;
1055     var1[5] = jet->EffectiveAreaCharged();
1056     var1[7] = ptLead;
1057     var1[8] = CheckAcceptance(phiJet,etaJet);
1058
1059     var2[1] = ptJet;
1060     var2[3] = etaJet;
1061     var2[4] = phiJet;
1062     var2[5] = jet->EffectiveAreaCharged();
1063     var2[7] = var1[8];
1064     var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1065
1066     if(ij<kMaxJets){
1067       fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1068       var1[0] = ij;
1069       var2[0] = ij;
1070       for(int it = 0;it <fNTrigger;it++){
1071         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1072           var1[6] = it;
1073           var2[6] = it;
1074           fhnJetPt[iType]->Fill(var1);
1075           fhnJetPtQA[iType]->Fill(var2);
1076         }
1077       }
1078     }
1079     var1[0] = kMaxJets;// fill for all jets
1080     var2[0] = kMaxJets;// fill for all jets
1081     for(int it = 0;it <fNTrigger;it++){
1082       if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1083         var1[6] = it;
1084         var2[6] = it;
1085         fhnJetPt[iType]->Fill(var1);
1086         fhnJetPtQA[iType]->Fill(var2);
1087       }
1088     }
1089
1090     fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
1091
1092     if(particlesList.GetSize()&&ij<kMaxJets){
1093       // Particles... correlated with jets...
1094       for(int it = 0;it<particlesList.GetEntries();++it){
1095         AliVParticle *part = (AliVParticle*)particlesList.At(it);
1096         Float_t deltaR = jet->DeltaR(part);
1097         if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1098       }
1099       // fill the jet shapes
1100     }// if we have particles
1101   }// Jet Loop
1102
1103
1104   // do something with dijets...
1105   if(ij0>=0&&ij1>0){
1106     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1107     Double_t ptJet0 = jet0->Pt();
1108     Double_t phiJet0 = jet0->Phi();
1109     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1110
1111     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1112     Double_t ptJet1 = jet1->Pt();
1113     Double_t phiJet1 = jet1->Phi();
1114     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1115
1116     Float_t deltaPhi = phiJet0 - phiJet1;
1117     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1118     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1119     deltaPhi = TMath::Abs(deltaPhi);
1120     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1121
1122     Float_t asym = 9999;
1123     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1124     fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1125     fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1126     fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1127     Float_t minv = 2.*(jet0->P()*jet1->P()-
1128                        jet0->Px()*jet1->Px()-
1129                        jet0->Py()*jet1->Py()-
1130                        jet0->Pz()*jet1->Pz());    // assume mass == 0;
1131     if(minv<0)minv=0; // prevent numerical instabilities
1132     minv = TMath::Sqrt(minv);
1133     fh1DijetMinv[iType]->Fill(minv);
1134   }
1135
1136
1137
1138   // count down the jets above thrueshold
1139   Int_t nOver = nJets;
1140   if(nOver>0){
1141     TIterator *jetIter = jetsList.MakeIterator();
1142     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1143     if(tmpJet){
1144       Float_t pt = tmpJet->Pt();
1145       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1146         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1147         while(pt<ptCut&&tmpJet){
1148           nOver--;
1149           tmpJet = (AliAODJet*)(jetIter->Next());
1150           if(tmpJet){
1151             pt = tmpJet->Pt();
1152           }
1153         }
1154         if(nOver<=0)break;
1155         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1156       }
1157     }
1158     delete jetIter;
1159   }
1160 }
1161
1162 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1163
1164   if(fFlagJetType[iType]<=0)return;
1165   Int_t refMult = fMultRec;
1166   if(iType==kJetGen||iType==kJetGenFull){
1167     refMult = fMultGen;
1168
1169   }
1170
1171   //
1172   Double_t var3[7]; // track number;p_{T};cent;#tracks;RP;cahrge
1173   var3[2] = fCentrality;
1174   var3[3] = refMult;
1175   Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
1176   var4[2] = fCentrality;
1177   Int_t nTrackOver = particlesList.GetSize();
1178   // do the same for tracks and jets
1179   if(nTrackOver>0){
1180     TIterator *trackIter = particlesList.MakeIterator();
1181     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1182     Float_t pt = tmpTrack->Pt();
1183     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1184       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1185       while(pt<ptCut&&tmpTrack){
1186         nTrackOver--;
1187         tmpTrack = (AliVParticle*)(trackIter->Next());
1188         if(tmpTrack){
1189           pt = tmpTrack->Pt();
1190         }
1191       }
1192       if(nTrackOver<=0)break;
1193       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1194     }
1195
1196
1197     trackIter->Reset();
1198     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1199     Float_t sumPt = 0;
1200
1201     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1202       Float_t tmpPt = tmpTrack->Pt();
1203       fh1PtTracksIn[iType]->Fill(tmpPt);
1204       fh1PtTracksInLow[iType]->Fill(tmpPt);
1205
1206       sumPt += tmpPt;
1207       Float_t tmpPhi = tmpTrack->Phi();
1208       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1209
1210
1211       Float_t phiRP = tmpPhi-fRPAngle;
1212       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1213       if(phiRP<0)phiRP += TMath::Pi();
1214       if(phiRP<0)phiRP += TMath::Pi();
1215       const float allPhi = -1./180.*TMath::Pi();
1216
1217       if(tmpPt<100){
1218         fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1219         fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1220
1221         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1222         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1223       }
1224       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1225       var3[0] = 1;
1226       var3[1] = tmpPt;
1227       var3[4] = phiBin;
1228       var3[6] = tmpTrack->Charge();
1229
1230       var4[0] = 1;
1231       var4[1] = tmpPt;
1232       var4[3] = tmpTrack->Eta();
1233       var4[4] = tmpPhi;
1234       var4[5] = tmpTrack->Charge();
1235
1236
1237       for(int it = 0;it <fNTrigger;it++){
1238         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1239           var3[0] = 1;
1240           var3[5] = it;
1241           fhnTrackPt[iType]->Fill(var3);
1242           if(tmpTrack==leading){
1243             var3[0] = 0;
1244             fhnTrackPt[iType]->Fill(var3);
1245           }
1246         }
1247       }
1248
1249
1250
1251
1252       fhnTrackPtQA[iType]->Fill(var4);
1253       if(tmpTrack==leading){
1254         var4[0] = 0;
1255         fhnTrackPtQA[iType]->Fill(var4);
1256         continue;
1257       }
1258
1259
1260
1261
1262
1263     }
1264     fh1SumPtTrack[iType]->Fill(sumPt);
1265     delete trackIter;
1266   }
1267
1268 }
1269
1270
1271 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1272
1273
1274   // Fill al the matching histos
1275   // It is important that the acceptances for the mathing are as large as possible
1276   // to avoid false matches on the edge of acceptance
1277   // therefore we add some extra matching jets as overhead
1278
1279   static TArrayI aGenIndex(recJetsList.GetEntries());
1280   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1281
1282   static TArrayI aRecIndex(genJetsList.GetEntries());
1283   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1284
1285   if(fDebug){
1286     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1287     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1288   }
1289   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1290                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1291                                             aGenIndex,aRecIndex,fDebug);
1292
1293   if(fDebug){
1294     for(int i = 0;i< aGenIndex.GetSize();++i){
1295       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1296     }
1297     for(int i = 0;i< aRecIndex.GetSize();++i){
1298       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1299     }
1300   }
1301
1302   Double_t container[8];
1303   Double_t containerRec[4];
1304   Double_t containerGen[4];
1305
1306   // loop over generated jets
1307   // consider the
1308   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1309     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1310     Double_t ptGen = genJet->Pt();
1311     Double_t phiGen = genJet->Phi();
1312     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1313     Double_t etaGen = genJet->Eta();
1314     containerGen[0] = ptGen;
1315     containerGen[1] = etaGen;
1316     containerGen[2] = phiGen;
1317     containerGen[3] = genJet->GetPtLeading();
1318
1319     fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1320
1321     if(JetSelected(genJet))
1322       fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1323
1324     Int_t ir = aRecIndex[ig];
1325     if(ir>=0&&ir<recJetsList.GetEntries()){
1326       AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1327
1328       fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1329
1330       if(JetSelected(genJet)){
1331         fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1332
1333         if(JetSelected(recJet)) {
1334           fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1335
1336         }
1337         containerGen[3] = recJet->GetPtLeading();
1338         fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner with leading track on reconstructed level
1339       }
1340     }
1341   }// loop over generated jets used for matching...
1342
1343
1344
1345   // loop over reconstructed jets
1346   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1347     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1348     Double_t etaRec = recJet->Eta();
1349     Double_t ptRec = recJet->Pt();
1350     Double_t phiRec = recJet->Phi();
1351     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1352
1353     containerRec[0] = ptRec;
1354     containerRec[1] = etaRec;
1355     containerRec[2] = phiRec;
1356     containerRec[3] = recJet->GetPtLeading();
1357
1358     container[0] = ptRec;
1359     container[1] = etaRec;
1360     container[2] = phiRec;
1361     container[3] = recJet->GetPtLeading();
1362
1363     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1364
1365     if(JetSelected(recJet)){
1366       fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
1367       // Fill Correlation
1368       Int_t ig = aGenIndex[ir];
1369       if(ig>=0 && ig<genJetsList.GetEntries()){
1370         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1371         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1372         Double_t ptGen  = genJet->Pt();
1373         Double_t phiGen = genJet->Phi();
1374         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1375         Double_t etaGen = genJet->Eta();
1376
1377         containerGen[0] = ptGen;
1378         containerGen[1] = etaGen;
1379         containerGen[2] = phiGen;
1380         containerGen[3] = genJet->GetPtLeading();
1381
1382         container[4] = ptGen;
1383         container[5] = etaGen;
1384         container[6] = phiGen;
1385         container[7] = genJet->GetPtLeading();
1386
1387         fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1388
1389         fhnCorrelation->Fill(container);
1390         if(ptGen>0){
1391           Float_t delta = (ptRec-ptGen)/ptGen;
1392           fh2RelPtFGen->Fill(ptGen,delta);
1393           fh3RelPtFGenLeadTrkPt->Fill(ptGen,delta,recJet->GetPtLeading());
1394           fh2PtFGen->Fill(ptGen,ptRec);
1395         }
1396       }
1397     }// loop over reconstructed jets
1398   }
1399   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1400 }
1401
1402
1403 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1404   //
1405   // Create the particle container for the correction framework manager and
1406   // link it
1407   //
1408   const Int_t kNvar   = 4 ; //number of variables on the grid:pt,eta, phi
1409   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1410   const Double_t kEtamin = -1.0, kEtamax = 1.0;
1411   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1412   const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
1413
1414   // can we neglect migration in eta and phi?
1415   // phi should be no problem since we cover full phi and are phi symmetric
1416   // eta migration is more difficult  due to needed acceptance correction
1417   // in limited eta range
1418
1419   //arrays for the number of bins in each dimension
1420   Int_t iBin[kNvar];
1421   iBin[0] = 125; //bins in pt
1422   iBin[1] =  4; //bins in eta
1423   iBin[2] =  1; // bins in phi
1424   iBin[3] =  10; // bins in leading track Pt
1425
1426   //arrays for lower bounds :
1427   Double_t* binEdges[kNvar];
1428   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1429     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1430
1431   //values for bin lower bounds
1432   //  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);
1433   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1434   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1435   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1436   binEdges[3][0]= kPtLeadingTrackPtMin;
1437   binEdges[3][1]= 1.;
1438   binEdges[3][2]= 2.;
1439   binEdges[3][3]= 3.;
1440   binEdges[3][4]= 4.;
1441   binEdges[3][5]= 5.;
1442   binEdges[3][6]= 6.;
1443   binEdges[3][7]= 8.;
1444   binEdges[3][8]= 10.;
1445   binEdges[3][9]= 12.;
1446   binEdges[3][10]= kPtLeadingTrackPtMax;
1447
1448   fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
1449   for (int k=0; k<kNvar; k++) {
1450     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1451   }
1452
1453   //create correlation matrix for unfolding
1454   Int_t thnDim[2*kNvar];
1455   for (int k=0; k<kNvar; k++) {
1456     //first half  : reconstructed
1457     //second half : MC
1458     thnDim[k]      = iBin[k];
1459     thnDim[k+kNvar] = iBin[k];
1460   }
1461
1462   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
1463   for (int k=0; k<kNvar; k++) {
1464     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1465     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1466   }
1467
1468   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1469     delete [] binEdges[ivar];
1470
1471
1472 }
1473
1474 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1475 {
1476   // Terminate analysis
1477   //
1478   if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1479 }
1480
1481
1482 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1483
1484   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1485
1486   Int_t iCount = 0;
1487   if(type==kTrackAOD){
1488     AliAODEvent *aod = 0;
1489     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1490     else aod = AODEvent();
1491     if(!aod){
1492       return iCount;
1493     }
1494     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1495       AliAODTrack *tr = aod->GetTrack(it);
1496       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1497       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1498       if(tr->Pt()<fMinTrackPt)continue;
1499       if(fDebug>0){
1500         if(tr->Pt()>20){
1501           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1502           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1503           tr->Print();
1504           //    tr->Dump();
1505           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1506           if(fESD){
1507             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1508             esdTr->Print("");
1509             // esdTr->Dump();
1510           }
1511         }
1512       }
1513       list->Add(tr);
1514       iCount++;
1515     }
1516   }
1517   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1518     AliMCEvent* mcEvent = MCEvent();
1519     if(!mcEvent)return iCount;
1520     // we want to have alivpartilces so use get track
1521     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1522       if(!mcEvent->IsPhysicalPrimary(it))continue;
1523       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1524       if(part->Pt()<fMinTrackPt)continue;
1525       if(type == kTrackKineAll){
1526         list->Add(part);
1527         iCount++;
1528       }
1529       else if(type == kTrackKineCharged){
1530         if(part->Particle()->GetPDG()->Charge()==0)continue;
1531         list->Add(part);
1532         iCount++;
1533       }
1534     }
1535   }
1536   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1537     AliAODEvent *aod = 0;
1538     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1539     else aod = AODEvent();
1540     if(!aod)return iCount;
1541     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1542     if(!tca)return iCount;
1543     for(int it = 0;it < tca->GetEntriesFast();++it){
1544       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1545       if(!part)continue;
1546       if(part->Pt()<fMinTrackPt)continue;
1547       if(!part->IsPhysicalPrimary())continue;
1548       if(type == kTrackAODMCAll){
1549         list->Add(part);
1550         iCount++;
1551       }
1552       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1553         if(part->Charge()==0)continue;
1554         if(kTrackAODMCCharged){
1555           list->Add(part);
1556         }
1557         else {
1558           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1559           list->Add(part);
1560         }
1561         iCount++;
1562       }
1563     }
1564   }// AODMCparticle
1565   list->Sort();
1566   return iCount;
1567
1568 }
1569
1570
1571 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1572   AliAODEvent *aod = 0;
1573   if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1574   else aod = AODEvent();
1575   if(!aod){
1576     return 101;
1577   }
1578   return aod->GetHeader()->GetCentrality();
1579 }
1580
1581
1582
1583 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1584   Bool_t selected = false;
1585
1586   if(!jet)return selected;
1587
1588   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1589     selected = kTRUE;
1590   }
1591   return selected;
1592
1593 }
1594
1595 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1596
1597   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1598   Int_t iCount = 0;
1599
1600   if(!jarray){
1601     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1602     return iCount;
1603   }
1604
1605
1606   for(int ij=0;ij<jarray->GetEntries();ij++){
1607     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1608     if(!jet)continue;
1609     if(type==0){
1610       // No cut at all, main purpose here is sorting
1611       list->Add(jet);
1612       iCount++;
1613     }
1614     else if(type == 1){
1615       // eta cut
1616       if(JetSelected(jet)){
1617         list->Add(jet);
1618         iCount++;
1619       }
1620     }
1621   }
1622
1623   list->Sort();
1624   return iCount;
1625
1626 }
1627
1628
1629 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1630   if(!jets)return 0;
1631
1632   Int_t refMult = 0;
1633   for(int ij = 0;ij < jets->GetEntries();++ij){
1634     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1635     if(!jet)continue;
1636     TRefArray *refs = jet->GetRefTracks();
1637     if(!refs)continue;
1638     refMult += refs->GetEntries();
1639   }
1640   return refMult;
1641
1642 }
1643
1644
1645 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1646   if(!jet)return 0;
1647   TRefArray *refs = jet->GetRefTracks();
1648   if(!refs) return 0;
1649   AliVParticle *leading = 0;
1650   Float_t fMaxPt = 0;
1651   for(int i = 0;i<refs->GetEntries();i++){
1652     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1653     if(!tmp)continue;
1654     if(tmp->Pt()>fMaxPt){
1655       leading = tmp;
1656       fMaxPt = tmp->Pt();
1657     }
1658   }
1659   return leading;
1660 }
1661
1662
1663 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1664   if(!jet)return 0;
1665   if(!list) return 0;
1666   AliVParticle *leading = 0;
1667   Float_t fMaxPt = 0;
1668   for(int i = 0;i<list->GetEntries();i++){
1669     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1670     if(!tmp)continue;
1671     if(jet->DeltaR(tmp)>r)continue;
1672     if(tmp->Pt()>fMaxPt){
1673       leading = tmp;
1674       fMaxPt = tmp->Pt();
1675     }
1676   }
1677   return leading;
1678 }
1679
1680 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1681 {
1682   Int_t phibin=-1;
1683   if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1684   Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1685   phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1686   if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1687   return phibin;
1688 }
1689
1690 void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1691   //
1692   // set number of trigger bins
1693   //
1694   if(n>0){
1695     fNTrigger = n;
1696     delete [] fTriggerName;
1697     fTriggerName = new TString [fNTrigger];
1698     delete [] fTriggerBit;fTriggerBit = 0;
1699     fTriggerBit = new UInt_t [fNTrigger];
1700   }
1701   else{
1702     fNTrigger = 0;
1703   }
1704 }
1705
1706
1707 void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1708   //
1709   // set trigger bin
1710   //
1711   if(i<fNTrigger){
1712     fTriggerBit[i] = it;
1713     fTriggerName[i] =  c;
1714   }
1715 }
1716
1717
1718 void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1719   //
1720   // Set number of acceptance bins
1721   //
1722
1723   if(n>0){
1724     fNAcceptance = n;
1725     delete [] fAcceptancePhiMin;
1726     delete [] fAcceptancePhiMax;
1727     delete [] fAcceptanceEtaMin;
1728     delete [] fAcceptanceEtaMax;
1729
1730     fAcceptancePhiMin = new Float_t[fNAcceptance];
1731     fAcceptancePhiMax = new Float_t[fNAcceptance];
1732     fAcceptanceEtaMin = new Float_t[fNAcceptance];
1733     fAcceptanceEtaMax = new Float_t[fNAcceptance];
1734   }
1735   else{
1736     fNTrigger = 0;
1737   }
1738 }
1739
1740 void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1741   //
1742   // Set acceptance bins
1743   //
1744
1745   if(i<fNAcceptance){
1746     Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1747
1748     fAcceptancePhiMin[i] = phiMin;
1749     fAcceptancePhiMax[i] = phiMax;
1750     fAcceptanceEtaMin[i] = etaMin;
1751     fAcceptanceEtaMax[i] = etaMax;
1752   }
1753   else{
1754     fNTrigger = 0;
1755   }
1756 }
1757
1758
1759 Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1760
1761   //
1762   // return aceptnace bin do no use ovelapping bins
1763   //
1764
1765   for(int i = 0;i<fNAcceptance;i++){
1766     if(phi<fAcceptancePhiMin[i])continue;
1767     if(phi>fAcceptancePhiMax[i])continue;
1768     if(eta<fAcceptanceEtaMin[i])continue;
1769     if(eta>fAcceptanceEtaMax[i])continue;
1770     return i;
1771   }
1772   // catch the rest
1773   return fNAcceptance;
1774 }
1775
1776 Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){
1777
1778   // invert the correction
1779   AliAODJet *jet = (AliAODJet*)list.At(0);  // highest pt jet
1780   if(!jet)return -1;
1781   if(jet->EffectiveAreaCharged()<=0)return -1;
1782   Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged();
1783   return rho;
1784
1785 }
1786
1787
1788
1789 AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1790   //
1791   delete [] fTriggerBit;
1792   delete [] fTriggerName;
1793
1794   delete [] fAcceptancePhiMin;
1795   delete [] fAcceptancePhiMax;
1796   delete [] fAcceptanceEtaMin;
1797   delete [] fAcceptanceEtaMax;
1798
1799
1800
1801 }