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