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