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 = 8;
 535     Int_t nBins2[nBinsSparse2] = {     kMaxJets+1, 25,   8,  18,             180, 10,fNTrigger,fNAcceptance+0.5};
 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};
 540     const Double_t xmax2[nBinsSparse2]  = {kMaxJets+0.5,250, 80, 0.9, 2.*TMath::Pi(),1.0,fNTrigger-0.5,fNAcceptance+0.5};
 541     fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin",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 = 5;
 570       const Int_t nBins4[nBinsSparse4] =    {    2, 50,  10,  20, 180};
 571       const Double_t xmin4[nBinsSparse4]  = { -0.5,  0,   0, -1.0,   0.};
 572       const Double_t xmax4[nBinsSparse4]  = {  1.5,150, 100,  1.0,2.*TMath::Pi()};
 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",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[8] = {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       if(ij<kMaxJets){
1030         if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
1031         var1[0] = ij;
1032         var2[0] = ij;
1033         for(int it = 0;it <fNTrigger;it++){
1034           if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1035             var1[6] = it;
1036             var2[6] = it;
1037             fhnJetPt[iType]->Fill(var1);
1038             fhnJetPtQA[iType]->Fill(var2);
1039           }
1040         }
1041       }
1042       var1[0] = kMaxJets;// fill for all jets
1043       var2[0] = kMaxJets;// fill for all jets
1044       for(int it = 0;it <fNTrigger;it++){
1045         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1046           var1[6] = it;
1047           fhnJetPt[iType]->Fill(var1);
1048           fhnJetPtQA[iType]->Fill(var2);
1049         }
1050       }
1051
1052       if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
1053
1054       if(particlesList.GetSize()&&ij<kMaxJets){
1055         // Particles... correlated with jets...
1056         for(int it = 0;it<particlesList.GetEntries();++it){
1057           AliVParticle *part = (AliVParticle*)particlesList.At(it);
1058           Float_t deltaR = jet->DeltaR(part);
1059           if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1060         }
1061         // fill the jet shapes
1062       }// if we have particles
1063   }// Jet Loop
1064
1065
1066   // do something with dijets...
1067   if(ij0>=0&&ij1>0){
1068     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1069     Double_t ptJet0 = jet0->Pt();
1070     Double_t phiJet0 = jet0->Phi();
1071     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1072
1073     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1074     Double_t ptJet1 = jet1->Pt();
1075     Double_t phiJet1 = jet1->Phi();
1076     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1077
1078     Float_t deltaPhi = phiJet0 - phiJet1;
1079     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1080     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1081     deltaPhi = TMath::Abs(deltaPhi);
1082     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1083
1084     Float_t asym = 9999;
1085     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1086       fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1087       fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1088       fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1089       Float_t minv = 2.*(jet0->P()*jet1->P()-
1090                          jet0->Px()*jet1->Px()-
1091                          jet0->Py()*jet1->Py()-
1092                          jet0->Pz()*jet1->Pz());    // assume mass == 0;
1093       if(minv<0)minv=0; // prevent numerical instabilities
1094       minv = TMath::Sqrt(minv);
1095       fh1DijetMinv[iType]->Fill(minv);
1096   }
1097
1098
1099
1100   // count down the jets above thrueshold
1101   Int_t nOver = nJets;
1102   if(nOver>0){
1103     TIterator *jetIter = jetsList.MakeIterator();
1104     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1105     if(tmpJet){
1106       Float_t pt = tmpJet->Pt();
1107       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1108         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1109         while(pt<ptCut&&tmpJet){
1110           nOver--;
1111           tmpJet = (AliAODJet*)(jetIter->Next());
1112           if(tmpJet){
1113             pt = tmpJet->Pt();
1114           }
1115         }
1116         if(nOver<=0)break;
1117         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1118       }
1119     }
1120     delete jetIter;
1121   }
1122 }
1123
1124 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1125
1126   if(fFlagJetType[iType]<=0)return;
1127   Int_t refMult = fMultRec;
1128   if(iType==kJetGen||iType==kJetGenFull){
1129     refMult = fMultGen;
1130
1131   }
1132
1133   //
1134   Double_t var3[6]; // track number;p_{T};cent;#tracks;RP
1135   var3[2] = fCentrality;
1136   var3[3] = refMult;
1137   Double_t var4[5]; // track number;p_{T};cent;#eta;#phi
1138   var4[2] = fCentrality;
1139   Int_t nTrackOver = particlesList.GetSize();
1140   // do the same for tracks and jets
1141   if(nTrackOver>0){
1142     TIterator *trackIter = particlesList.MakeIterator();
1143     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1144     Float_t pt = tmpTrack->Pt();
1145     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1146       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1147       while(pt<ptCut&&tmpTrack){
1148         nTrackOver--;
1149         tmpTrack = (AliVParticle*)(trackIter->Next());
1150         if(tmpTrack){
1151           pt = tmpTrack->Pt();
1152         }
1153       }
1154       if(nTrackOver<=0)break;
1155       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1156     }
1157
1158
1159     trackIter->Reset();
1160     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1161     Float_t sumPt = 0;
1162
1163     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1164       Float_t tmpPt = tmpTrack->Pt();
1165       fh1PtTracksIn[iType]->Fill(tmpPt);
1166       fh1PtTracksInLow[iType]->Fill(tmpPt);
1167
1168       sumPt += tmpPt;
1169       Float_t tmpPhi = tmpTrack->Phi();
1170       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1171
1172
1173       Float_t phiRP = tmpPhi-fRPAngle;
1174       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1175       if(phiRP<0)phiRP += TMath::Pi();
1176       if(phiRP<0)phiRP += TMath::Pi();
1177       const float allPhi = -1./180.*TMath::Pi();
1178
1179       if(tmpPt<100){
1180         fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1181         fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1182
1183         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1184         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1185       }
1186       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1187       var3[0] = 1;
1188       var3[1] = tmpPt;
1189       var3[4] = phiBin;
1190
1191       var4[0] = 1;
1192       var4[1] = tmpPt;
1193       var4[3] = tmpTrack->Eta();
1194       var4[4] = tmpPhi;
1195
1196
1197
1198       for(int it = 0;it <fNTrigger;it++){
1199         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1200           var3[0] = 1;
1201           var3[5] = it;
1202           fhnTrackPt[iType]->Fill(var3);
1203           if(tmpTrack==leading){
1204             var3[0] = 0;
1205             fhnTrackPt[iType]->Fill(var3);
1206           }
1207         }
1208       }
1209
1210
1211
1212
1213
1214       fhnTrackPtQA[iType]->Fill(var4);
1215
1216       if(tmpTrack==leading){
1217         var4[0] = 0;
1218         fhnTrackPtQA[iType]->Fill(var4);
1219         continue;
1220       }
1221     }
1222     fh1SumPtTrack[iType]->Fill(sumPt);
1223     delete trackIter;
1224   }
1225
1226 }
1227
1228
1229 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1230
1231
1232   // Fill al the matching histos
1233   // It is important that the acceptances for the mathing are as large as possible
1234   // to avoid false matches on the edge of acceptance
1235   // therefore we add some extra matching jets as overhead
1236
1237   static TArrayI aGenIndex(recJetsList.GetEntries());
1238   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1239
1240   static TArrayI aRecIndex(genJetsList.GetEntries());
1241   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1242
1243   if(fDebug){
1244     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1245     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1246   }
1247   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1248                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1249                                             aGenIndex,aRecIndex,fDebug);
1250
1251   if(fDebug){
1252     for(int i = 0;i< aGenIndex.GetSize();++i){
1253       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1254     }
1255     for(int i = 0;i< aRecIndex.GetSize();++i){
1256       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1257     }
1258   }
1259
1260   Double_t container[6];
1261
1262   // loop over generated jets
1263   // consider the
1264   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1265     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1266     Double_t ptGen = genJet->Pt();
1267     Double_t phiGen = genJet->Phi();
1268     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1269     Double_t etaGen = genJet->Eta();
1270     container[3] = ptGen;
1271     container[4] = etaGen;
1272     container[5] = phiGen;
1273     fhnJetContainer->Fill(&container[3],kStep0);
1274     if(JetSelected(genJet)){
1275       fhnJetContainer->Fill(&container[3],kStep1);
1276       Int_t ir = aRecIndex[ig];
1277       if(ir>=0&&ir<recJetsList.GetEntries()){
1278         fhnJetContainer->Fill(&container[3],kStep2);
1279         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1280         if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep4);
1281         if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep3);
1282       }
1283     }
1284   }// loop over generated jets used for matching...
1285
1286
1287
1288   // loop over reconstructed jets
1289   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1290     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1291     Double_t etaRec = recJet->Eta();
1292     Double_t ptRec = recJet->Pt();
1293     Double_t phiRec = recJet->Phi();
1294     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1295     // do something with dijets...
1296
1297     container[0] = ptRec;
1298     container[1] = etaRec;
1299     container[2] = phiRec;
1300
1301     fhnJetContainer->Fill(container,kStep0+kMaxStep);
1302     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1303
1304     if(JetSelected(recJet)){
1305       fhnJetContainer->Fill(container,kStep1+kMaxStep);
1306       // Fill Correlation
1307       Int_t ig = aGenIndex[ir];
1308       if(ig>=0 && ig<genJetsList.GetEntries()){
1309         fhnJetContainer->Fill(container,kStep2+kMaxStep);
1310         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1311         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1312         Double_t ptGen  = genJet->Pt();
1313         Double_t phiGen = genJet->Phi();
1314         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1315         Double_t etaGen = genJet->Eta();
1316
1317         container[3] = ptGen;
1318         container[4] = etaGen;
1319         container[5] = phiGen;
1320         //
1321         // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1322         //
1323         if(JetSelected(genJet))fhnJetContainer->Fill(container,kStep4+kMaxStep);
1324         fhnJetContainer->Fill(container,kStep3+kMaxStep);
1325         fhnCorrelation->Fill(container,0);
1326         if(ptGen>0){
1327           Float_t delta = (ptRec-ptGen)/ptGen;
1328           fh2RelPtFGen->Fill(ptGen,delta);
1329           fh2PtFGen->Fill(ptGen,ptRec);
1330         }
1331       }
1332     }// loop over reconstructed jets
1333   }
1334   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1335 }
1336
1337
1338 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1339   //
1340   // Create the particle container for the correction framework manager and
1341   // link it
1342   //
1343   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1344   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1345   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1346   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1347
1348   // can we neglect migration in eta and phi?
1349   // phi should be no problem since we cover full phi and are phi symmetric
1350   // eta migration is more difficult  due to needed acceptance correction
1351   // in limited eta range
1352
1353   //arrays for the number of bins in each dimension
1354   Int_t iBin[kNvar];
1355   iBin[0] = 125; //bins in pt
1356   iBin[1] =  1; //bins in eta
1357   iBin[2] = 1; // bins in phi
1358
1359   //arrays for lower bounds :
1360   Double_t* binEdges[kNvar];
1361   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1362     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1363
1364   //values for bin lower bounds
1365   //  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);
1366   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1367   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1368   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1369
1370
1371   fhnJetContainer = new AliTHn(Form("fahnJetContainer"),Form("AliTHn jet info"),kMaxStep*2,kNvar,iBin);
1372   for (int k=0; k<kNvar; k++) {
1373     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1374   }
1375
1376   //create correlation matrix for unfolding
1377   Int_t thnDim[2*kNvar];
1378   for (int k=0; k<kNvar; k++) {
1379     //first half  : reconstructed
1380     //second half : MC
1381     thnDim[k]      = iBin[k];
1382     thnDim[k+kNvar] = iBin[k];
1383   }
1384
1385   fhnCorrelation = new AliTHn("fahnCorrelation","AliTHn with correlations",1,2*kNvar,thnDim);
1386   for (int k=0; k<kNvar; k++) {
1387     fhnCorrelation->SetBinLimits(k,binEdges[k]);
1388     fhnCorrelation->SetBinLimits(k+kNvar,binEdges[k]);
1389   }
1390
1391   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1392     delete [] binEdges[ivar];
1393
1394
1395 }
1396
1397 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1398 {
1399 // Terminate analysis
1400 //
1401     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1402 }
1403
1404
1405 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1406
1407   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1408
1409   Int_t iCount = 0;
1410   if(type==kTrackAOD){
1411     AliAODEvent *aod = 0;
1412     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1413     else aod = AODEvent();
1414     if(!aod){
1415       return iCount;
1416     }
1417     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1418       AliAODTrack *tr = aod->GetTrack(it);
1419       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1420       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1421       if(tr->Pt()<fMinTrackPt)continue;
1422       if(fDebug>0){
1423         if(tr->Pt()>20){
1424           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1425           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1426           tr->Print();
1427           //    tr->Dump();
1428           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1429           if(fESD){
1430             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1431             esdTr->Print("");
1432             // esdTr->Dump();
1433           }
1434         }
1435       }
1436       list->Add(tr);
1437       iCount++;
1438     }
1439   }
1440   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1441     AliMCEvent* mcEvent = MCEvent();
1442     if(!mcEvent)return iCount;
1443     // we want to have alivpartilces so use get track
1444     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1445       if(!mcEvent->IsPhysicalPrimary(it))continue;
1446       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1447       if(part->Pt()<fMinTrackPt)continue;
1448       if(type == kTrackKineAll){
1449         list->Add(part);
1450         iCount++;
1451       }
1452       else if(type == kTrackKineCharged){
1453         if(part->Particle()->GetPDG()->Charge()==0)continue;
1454         list->Add(part);
1455         iCount++;
1456       }
1457     }
1458   }
1459   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1460     AliAODEvent *aod = 0;
1461     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1462     else aod = AODEvent();
1463     if(!aod)return iCount;
1464     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1465     if(!tca)return iCount;
1466     for(int it = 0;it < tca->GetEntriesFast();++it){
1467       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1468       if(!part)continue;
1469       if(part->Pt()<fMinTrackPt)continue;
1470       if(!part->IsPhysicalPrimary())continue;
1471       if(type == kTrackAODMCAll){
1472         list->Add(part);
1473         iCount++;
1474       }
1475       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1476         if(part->Charge()==0)continue;
1477         if(kTrackAODMCCharged){
1478           list->Add(part);
1479         }
1480         else {
1481           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1482           list->Add(part);
1483         }
1484         iCount++;
1485       }
1486     }
1487   }// AODMCparticle
1488   list->Sort();
1489   return iCount;
1490
1491 }
1492
1493
1494 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1495     AliAODEvent *aod = 0;
1496     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1497     else aod = AODEvent();
1498     if(!aod){
1499       return 101;
1500     }
1501     return aod->GetHeader()->GetCentrality();
1502 }
1503
1504
1505
1506 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1507   Bool_t selected = false;
1508
1509   if(!jet)return selected;
1510
1511   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1512     selected = kTRUE;
1513   }
1514   return selected;
1515
1516 }
1517
1518 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1519
1520   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1521   Int_t iCount = 0;
1522
1523   if(!jarray){
1524     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1525     return iCount;
1526   }
1527
1528
1529   for(int ij=0;ij<jarray->GetEntries();ij++){
1530     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1531     if(!jet)continue;
1532     if(type==0){
1533       // No cut at all, main purpose here is sorting
1534       list->Add(jet);
1535       iCount++;
1536     }
1537     else if(type == 1){
1538       // eta cut
1539       if(JetSelected(jet)){
1540         list->Add(jet);
1541         iCount++;
1542       }
1543     }
1544   }
1545
1546   list->Sort();
1547   return iCount;
1548
1549 }
1550
1551
1552 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1553   if(!jets)return 0;
1554
1555   Int_t refMult = 0;
1556   for(int ij = 0;ij < jets->GetEntries();++ij){
1557     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1558     if(!jet)continue;
1559     TRefArray *refs = jet->GetRefTracks();
1560     if(!refs)continue;
1561     refMult += refs->GetEntries();
1562   }
1563   return refMult;
1564
1565 }
1566
1567
1568 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1569   if(!jet)return 0;
1570   TRefArray *refs = jet->GetRefTracks();
1571   if(!refs) return 0;
1572   AliVParticle *leading = 0;
1573   Float_t fMaxPt = 0;
1574   for(int i = 0;i<refs->GetEntries();i++){
1575     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1576     if(!tmp)continue;
1577     if(tmp->Pt()>fMaxPt){
1578       leading = tmp;
1579       fMaxPt = tmp->Pt();
1580     }
1581   }
1582   return leading;
1583 }
1584
1585
1586 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1587   if(!jet)return 0;
1588   if(!list) return 0;
1589   AliVParticle *leading = 0;
1590   Float_t fMaxPt = 0;
1591   for(int i = 0;i<list->GetEntries();i++){
1592     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1593     if(!tmp)continue;
1594     if(jet->DeltaR(tmp)>r)continue;
1595     if(tmp->Pt()>fMaxPt){
1596       leading = tmp;
1597       fMaxPt = tmp->Pt();
1598     }
1599   }
1600   return leading;
1601 }
1602
1603 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1604 {
1605     Int_t phibin=-1;
1606     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1607     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1608     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1609     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1610     return phibin;
1611 }
1612
1613 void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1614   //
1615   // set number of trigger bins
1616   //
1617   if(n>0){
1618       fNTrigger = n;
1619       delete [] fTriggerName;
1620       fTriggerName = new TString [fNTrigger];
1621       delete [] fTriggerBit;fTriggerBit = 0;
1622       fTriggerBit = new UInt_t [fNTrigger];
1623   }
1624   else{
1625     fNTrigger = 0;
1626   }
1627 }
1628
1629
1630 void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1631   //
1632   // set trigger bin
1633   //
1634   if(i<fNTrigger){
1635     fTriggerBit[i] = it;
1636     fTriggerName[i] =  c;
1637   }
1638 }
1639
1640
1641 void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1642   //
1643   // Set number of acceptance bins
1644   //
1645
1646   if(n>0){
1647       fNAcceptance = n;
1648       delete [] fAcceptancePhiMin;
1649       delete [] fAcceptancePhiMax;
1650       delete [] fAcceptanceEtaMin;
1651       delete [] fAcceptanceEtaMax;
1652
1653       fAcceptancePhiMin = new Float_t[fNAcceptance];
1654       fAcceptancePhiMax = new Float_t[fNAcceptance];
1655       fAcceptanceEtaMin = new Float_t[fNAcceptance];
1656       fAcceptanceEtaMax = new Float_t[fNAcceptance];
1657   }
1658   else{
1659     fNTrigger = 0;
1660   }
1661 }
1662
1663 void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1664   //
1665   // Set acceptance bins
1666   //
1667
1668   if(i<fNAcceptance){
1669     Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1670
1671     fAcceptancePhiMin[i] = phiMin;
1672     fAcceptancePhiMax[i] = phiMax;
1673     fAcceptanceEtaMin[i] = etaMin;
1674     fAcceptanceEtaMax[i] = etaMax;
1675   }
1676   else{
1677     fNTrigger = 0;
1678   }
1679 }
1680
1681
1682 Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1683
1684   //
1685   // return aceptnace bin do no use ovelapping bins
1686   //
1687
1688   for(int i = 0;i<fNAcceptance;i++){
1689     if(phi<fAcceptancePhiMin[i])continue;
1690     if(phi>fAcceptancePhiMax[i])continue;
1691     if(eta<fAcceptanceEtaMin[i])continue;
1692     if(eta>fAcceptanceEtaMax[i])continue;
1693     return i;
1694   }
1695   // catch the rest
1696   return fNAcceptance;
1697 }
1698
1699
1700
1701
1702
1703 AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1704   //
1705   delete [] fTriggerBit;
1706   delete [] fTriggerName;
1707
1708   delete [] fAcceptancePhiMin;
1709   delete [] fAcceptancePhiMax;
1710   delete [] fAcceptanceEtaMin;
1711   delete [] fAcceptanceEtaMax;
1712
1713
1714
1715 }