PWGJE/UserTasks/AliAnalysisTaskJetHBOM.cxx

   1 // **************************************
   2 // Task used for the correction of detector effects for background fluctuations in jet spectra by the HBOM method
   3 // *******************************************
   4
   5
   6 /**************************************************************************
   7  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   8  *                                                                        *
   9  * Author: The ALICE Off-line Project.                                    *
  10  * Contributors are mentioned in the code where appropriate.              *
  11  *                                                                        *
  12  * Permission to use, copy, modify and distribute this software and its   *
  13  * documentation strictly for non-commercial purposes is hereby granted   *
  14  * without fee, provided that the above copyright notice appears in all   *
  15  * copies and that both the copyright notice and this permission notice   *
  16  * appear in the supporting documentation. The authors make no claims     *
  17  * about the suitability of this software for any purpose. It is          *
  18  * provided "as is" without express or implied warranty.                  *
  19  **************************************************************************/
  20
  21
  22 #include <TROOT.h>
  23 #include <TRandom3.h>
  24 #include <TSystem.h>
  25 #include <TInterpreter.h>
  26 #include <TChain.h>
  27 #include <TRefArray.h>
  28 #include <TFile.h>
  29 #include <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TF1.h>
  35 #include <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetHBOM.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDEvent.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODHandler.h"
  45 #include "AliAODExtension.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODJet.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliMCEvent.h"
  51 #include "AliStack.h"
  52 #include "AliGenPythiaEventHeader.h"
  53 #include "AliGenCocktailEventHeader.h"
  54 #include "AliInputEventHandler.h"
  55 #include "AliAODJetEventBackground.h"
  56
  57
  58 using std::cout;
  59 using std::endl;
  60 using std::vector;
  61
  62 ClassImp(AliAnalysisTaskJetHBOM)
  63
  64 AliAnalysisTaskJetHBOM::~AliAnalysisTaskJetHBOM(){
  65   //
  66   // Destructor
  67   //
  68
  69   delete fRef;
  70   delete fRandom;
  71
  72   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
  73   delete fTCARandomConesOut;
  74
  75   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
  76   delete fTCARandomConesOutRan;
  77
  78 }
  79
  80 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM():
  81   AliAnalysisTaskSE(),
  82   fAOD(0x0),
  83   fAODExtension(0x0),
  84   fRef(new TRefArray),
  85   fUseAODTrackInput(kFALSE),
  86   fUseAODMCInput(kFALSE),
  87   fEventSelection(kFALSE),
  88   fFilterMask(0),
  89   fFilterMaskBestPt(0),
  90   fFilterType(0),
  91   fJetTypes(kJet),
  92   fTrackTypeRec(kTrackUndef),
  93   fTrackTypeGen(kTrackUndef),
  94   fNSkipLeadingRan(0),
  95   fNSkipLeadingCone(0),
  96   fNRandomCones(0),
  97   randCone_pos(0),
  98   randCone_Eta(0),
  99   randCone_Phi(0),
 100   fNHBOM(0),
 101   fTrackEtaWindow(0.9),
 102   //fRecEtaWindow(0.5),
 103   fTrackPtCut(0.),
 104   fJetOutputMinPt(0.150),
 105   fMaxTrackPtInJet(100.),
 106   fVtxZCut(8),
 107   fVtxR2Cut(1),
 108   fCentCutUp(0),
 109   fCentCutLo(0),
 110   fNonStdBranch(""),
 111   fBackgroundBranch(""),
 112   fNonStdFile(""),
 113   fMomResH1(0x0),
 114   fMomResH2(0x0),
 115   fMomResH3(0x0),
 116   fMomResH1Fit(0x0),
 117   fMomResH2Fit(0x0),
 118   fMomResH3Fit(0x0),
 119   fhEffH1(0x0),
 120   fhEffH2(0x0),
 121   fhEffH3(0x0),
 122   fUseTrMomentumSmearing(kFALSE),
 123   fUseDiceEfficiency(kFALSE),
 124   fRparam(1.0),
 125   fAlgorithm(fastjet::kt_algorithm),
 126   fStrategy(fastjet::Best),
 127   fRecombScheme(fastjet::BIpt_scheme),
 128   fAreaType(fastjet::active_area),
 129   fGhostArea(0.01),
 130   fActiveAreaRepeats(1),
 131   fGhostEtamax(1.5),
 132   background(0),
 133   fTCARandomConesOut(0x0),
 134   fTCARandomConesOutRan(0x0),
 135   fRandom(0),
 136   fh1Xsec(0x0),
 137   fh1Trials(0x0),
 138   fh1PtHard(0x0),
 139   fh1PtHardNoW(0x0),
 140   fh1PtHardTrials(0x0),
 141   fh1Nch(0x0),
 142   fh1CentralityPhySel(0x0),
 143   fh1Centrality(0x0),
 144   fh1DeltapT(0x0),
 145   fh1Rho(0x0),
 146   fh1PtRandCone(0x0),
 147   fh1Area(0x0),
 148   fh1efficiencyPt(0x0),
 149   fh2efficiencyPhi(0x0),
 150   fh1ZPhySel(0x0),
 151   fh1Z(0x0),
 152   fHistList(0x0)
 153 {
 154
 155 }
 156
 157 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(const char* name):
 158   AliAnalysisTaskSE(name),
 159   fAOD(0x0),
 160   fAODExtension(0x0),
 161   fRef(new TRefArray),
 162   fUseAODTrackInput(kFALSE),
 163   fUseAODMCInput(kFALSE),
 164   fEventSelection(kFALSE),
 165   fFilterMask(0),
 166   fFilterMaskBestPt(0),
 167   fFilterType(0),
 168   fJetTypes(kJet),
 169   fTrackTypeRec(kTrackUndef),
 170   fTrackTypeGen(kTrackUndef),
 171   fNSkipLeadingRan(0),
 172   fNSkipLeadingCone(0),
 173   fNRandomCones(0),
 174   randCone_pos(0),
 175   randCone_Eta(0),
 176   randCone_Phi(0),
 177   fNHBOM(0),
 178   fTrackEtaWindow(0.9),
 179   //fRecEtaWindow(0.5),
 180   fTrackPtCut(0.),
 181   fJetOutputMinPt(0.150),
 182   fMaxTrackPtInJet(100.),
 183   fVtxZCut(8),
 184   fVtxR2Cut(1),
 185   fCentCutUp(0),
 186   fCentCutLo(0),
 187   fNonStdBranch(""),
 188   fBackgroundBranch(""),
 189   fNonStdFile(""),
 190   fMomResH1(0x0),
 191   fMomResH2(0x0),
 192   fMomResH3(0x0),
 193   fMomResH1Fit(0x0),
 194   fMomResH2Fit(0x0),
 195   fMomResH3Fit(0x0),
 196   fhEffH1(0x0),
 197   fhEffH2(0x0),
 198   fhEffH3(0x0),
 199   fUseTrMomentumSmearing(kFALSE),
 200   fUseDiceEfficiency(kFALSE),
 201   fRparam(1.0),
 202   fAlgorithm(fastjet::kt_algorithm),
 203   fStrategy(fastjet::Best),
 204   fRecombScheme(fastjet::BIpt_scheme),
 205   fAreaType(fastjet::active_area),
 206   fGhostArea(0.01),
 207   fActiveAreaRepeats(1),
 208   fGhostEtamax(1.5),
 209   background(0),
 210   fTCARandomConesOut(0x0),
 211   fTCARandomConesOutRan(0x0),
 212   fRandom(0),
 213   fh1Xsec(0x0),
 214   fh1Trials(0x0),
 215   fh1PtHard(0x0),
 216   fh1PtHardNoW(0x0),
 217   fh1PtHardTrials(0x0),
 218   fh1Nch(0x0),
 219   fh1CentralityPhySel(0x0),
 220   fh1Centrality(0x0),
 221   fh1DeltapT(0x0),
 222   fh1Rho(0x0),
 223   fh1PtRandCone(0x0),
 224   fh1Area(0x0),
 225   fh1efficiencyPt(0x0),
 226   fh2efficiencyPhi(0x0),
 227   fh1ZPhySel(0x0),
 228   fh1Z(0x0),
 229   fHistList(0x0)
 230 {
 231  DefineOutput(1, TList::Class());
 232 }
 233
 234
 235
 236 Bool_t AliAnalysisTaskJetHBOM::Notify()
 237 {
 238   //
 239   // Implemented Notify() to read the cross sections
 240   // and number of trials from pyxsec.root
 241   //
 242   return kTRUE;
 243 }
 244
 245 void AliAnalysisTaskJetHBOM::UserCreateOutputObjects()
 246 {
 247
 248   //
 249   // Create the output container
 250   //
 251
 252   fRandom = new TRandom3(0);
 253
 254
 255   // Connect the AOD
 256
 257
 258   if (fDebug > 1) printf("AnalysisTaskJetHBOM::UserCreateOutputObjects() \n");
 259
 260
 261
 262   if(fNonStdBranch.Length()!=0)
 263     {
 264       // only create the output branch if we have a name
 265       // Create a new branch for jets...
 266       //  -> cleared in the UserExec....
 267       // here we can also have the case that the brnaches are written to a separate file
 268
 269       // create the branch for the random cones with the same R
 270       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
 271       if(fUseDiceEfficiency || fUseTrMomentumSmearing)
 272         cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
 273
 274       //create array for the random cones; Until now only one cone per event is used
 275       if(!AODEvent()->FindListObject(cName.Data())){
 276         fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
 277         fTCARandomConesOut->SetName(cName.Data());
 278         AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
 279       }
 280       // create the branch with the random for the random cones on the random event
 281       cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 282       if(!AODEvent()->FindListObject(cName.Data())){
 283         fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
 284         fTCARandomConesOutRan->SetName(cName.Data());
 285         AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
 286       }
 287
 288       if(fNonStdFile.Length()!=0){
 289         //
 290         // case that we have an AOD extension we need to fetch the jets from the extended output
 291         // we identify the extension aod event by looking for the branchname
 292         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 293         // case that we have an AOD extension we need can fetch the background maybe from the extended output
 294         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 295       }
 296     }
 297
 298   //  FitMomentumResolution();
 299
 300
 301   if(!fHistList)fHistList = new TList();
 302   fHistList->SetOwner();
 303   PostData(1, fHistList); // post data in any case once
 304
 305   Bool_t oldStatus = TH1::AddDirectoryStatus();
 306   TH1::AddDirectory(kFALSE);
 307
 308   //
 309   //  Histogram
 310
 311   const Int_t nBinPt = 100;
 312   Double_t binLimitsPt[nBinPt+1];
 313   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 314     if(iPt == 0){
 315       binLimitsPt[iPt] = 0.0;
 316     }
 317     else {// 1.0
 318       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
 319     }
 320   }
 321
 322   const Int_t nBinPhi = 90;
 323   Double_t binLimitsPhi[nBinPhi+1];
 324   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 325     if(iPhi==0){
 326       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 327     }
 328     else{
 329       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 330     }
 331   }
 332
 333
 334
 335   const Int_t nBinEta = 40;
 336   Double_t binLimitsEta[nBinEta+1];
 337   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 338     if(iEta==0){
 339       binLimitsEta[iEta] = -2.0;
 340     }
 341     else{
 342       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 343     }
 344   }
 345
 346   const int nChMax = 5000;
 347
 348   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 349   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 350
 351   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 352   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 353
 354
 355   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 356   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 357   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 358
 359   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 360
 361   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 362   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 363
 364   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 365   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 366
 367   fh1DeltapT = new TH1F("fh1DeltapT","DeltapT",100,-50,50);
 368   fh1Rho = new TH1F("fh1Rho","Rho",100,0,200);
 369   fh1PtRandCone = new TH1F("fh1PtRandCone","pt",100,0,200);
 370   fh1Area = new TH1F("fh1Area","area",100,0,1);
 371
 372   const Int_t saveLevel = 3; // large save level more histos
 373   if(saveLevel>0){
 374     fHistList->Add(fh1Xsec);
 375     fHistList->Add(fh1Trials);
 376
 377     fHistList->Add(fh1Nch);
 378     fHistList->Add(fh1Centrality);
 379     fHistList->Add(fh1CentralityPhySel);
 380     fHistList->Add(fh1Z);
 381     fHistList->Add(fh1ZPhySel);
 382     fHistList->Add(fh1DeltapT);
 383     fHistList->Add(fh1Rho);
 384     fHistList->Add(fh1PtRandCone);
 385     fHistList->Add(fh1Area);
 386   }
 387
 388   // =========== Switch on Sumw2 for all histos ===========
 389   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 390     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 391     if (h1){
 392       h1->Sumw2();
 393       continue;
 394     }
 395     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 396     if(hn)hn->Sumw2();
 397   }
 398   TH1::AddDirectory(oldStatus);
 399 }
 400
 401 void AliAnalysisTaskJetHBOM::Init()
 402 {
 403   //
 404   // Initialization
 405   //
 406
 407   if (fDebug > 1) printf("AnalysisTaskJetHBOM::Init() \n");
 408
 409   FitMomentumResolution();
 410
 411 }
 412
 413 void AliAnalysisTaskJetHBOM::UserExec(Option_t */*option*/)
 414 {
 415
 416   // handle and reset the output jet branch
 417
 418   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
 419   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
 420
 421   //
 422   // Execute analysis for current event
 423   //
 424   AliESDEvent *fESD = 0;
 425   if(fUseAODTrackInput){
 426     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 427     if(!fAOD){
 428       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 429       return;
 430     }
 431     // fetch the header
 432   }
 433   else{
 434     //  assume that the AOD is in the general output...
 435     fAOD  = AODEvent();
 436     if(!fAOD){
 437       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 438       return;
 439     }
 440     if(fDebug>0){
 441       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 442     }
 443   }
 444
 445   //Check if information is provided detector level effects
 446   if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
 447   if(!fhEffH1 || !fhEffH2 || !fhEffH3)       fUseDiceEfficiency = kFALSE;
 448
 449   Bool_t selectEvent =  false;
 450   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 451
 452   Float_t cent = 0;
 453   Float_t zVtx  = 0;
 454
 455   if(fAOD){
 456     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 457     TString vtxTitle(vtxAOD->GetTitle());
 458     zVtx = vtxAOD->GetZ();
 459
 460     cent = fAOD->GetHeader()->GetCentrality();
 461     if(physicsSelection){
 462       fh1CentralityPhySel->Fill(cent);
 463       fh1ZPhySel->Fill(zVtx);
 464     }
 465     // zVertex and centrality selection
 466     if(fEventSelection){
 467       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 468         Float_t yvtx = vtxAOD->GetY();
 469         Float_t xvtx = vtxAOD->GetX();
 470         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 471         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){//usual fVtxZCut=10 and fVtxR2Cut=1  // apply vertex cut later on
 472           if(physicsSelection){
 473             selectEvent = true;
 474           }
 475         }
 476       }
 477       //centrality selection
 478       if(fCentCutUp>0){
 479         if(cent<fCentCutLo||cent>fCentCutUp){
 480           selectEvent = false;
 481         }
 482       }
 483     }else{
 484       selectEvent = true;
 485     }
 486   }
 487
 488
 489   if(!selectEvent){
 490     PostData(1, fHistList);
 491     return;
 492   }
 493   fh1Centrality->Fill(cent);
 494   fh1Z->Fill(zVtx);
 495   fh1Trials->Fill("#sum{ntrials}",1);
 496
 497
 498   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 499
 500   // ==== General variables needed
 501
 502
 503
 504   // we simply fetch the tracks/mc particles as a list of AliVParticles
 505
 506   //reconstructed particles
 507   TList recParticles;
 508   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 509   Float_t nCh = recParticles.GetEntries();
 510   fh1Nch->Fill(nCh);
 511   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 512   //nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 513   //if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 514
 515
 516   //apply efficiency fNHBOM times
 517   if(fNHBOM>0){
 518     for(int particle=0;particle<recParticles.GetEntries();particle++){
 519       // hier Effizienzen laden und überprüfen ob das Teilchen nachgewiesen wird.
 520       AliVParticle *vp = (AliVParticle*)recParticles.At(particle);
 521
 522       Double_t pT = vp->Pt();
 523       Double_t phi = vp->Phi();
 524
 525       //load efficiency
 526       Double_t efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
 527       Double_t efficiencyPhi = fh2efficiencyPhi->GetBinContent(fh2efficiencyPhi->FindBin(phi,pT));
 528       if(pT>10){
 529         efficiencyPt = 0.857; //this is the result for the fit with pT>10; statistic is low for pT>10
 530         //if the efficiency is not from fastMCInput_LHC10h_110719a.root this should be calculated new
 531         if(fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))>0.9||fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))<0.8){
 532           efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
 533           Printf("%s:%d Wrong efficiency input. Check efficiency for pt>10GeV",(char*)__FILE__,__LINE__);
 534         }
 535       }
 536       Double_t eff = efficiencyPt*efficiencyPhi; //over all efficiency
 537
 538       // if ran<eff -> particle is detected eff^fNHBOM = efficiency to detect particle fNHBOM times
 539       Double_t ran = fRandom->Rndm();
 540       if(ran>TMath::Power(eff,fNHBOM)){
 541         recParticles.Remove(vp);
 542       }
 543     }
 544   }
 545
 546
 547   // find the jets....
 548
 549   vector<fastjet::PseudoJet> inputParticlesRec;
 550   vector<fastjet::PseudoJet> inputParticlesRecRan;
 551
 552   //randomize particles
 553   AliAODJet vTmpRan(1,0,0,1);
 554   for(int i = 0; i < recParticles.GetEntries(); i++){
 555     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 556
 557     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 558     // we take total momentum here
 559
 560       //Add particles to fastjet in case we are not running toy model
 561       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 562       jInp.set_user_index(i);
 563       inputParticlesRec.push_back(jInp);
 564
 565     // the randomized input changes eta and phi, but keeps the p_T
 566       Double_t pT = vp->Pt();
 567       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
 568       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
 569
 570       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
 571       Double_t pZ = pT/TMath::Tan(theta);
 572
 573       Double_t pX = pT * TMath::Cos(phi);
 574       Double_t pY = pT * TMath::Sin(phi);
 575       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 576       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 577
 578       jInpRan.set_user_index(i);
 579       inputParticlesRecRan.push_back(jInpRan);
 580       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
 581
 582   }// recparticles
 583
 584   if(inputParticlesRec.size()==0){
 585     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 586     PostData(1, fHistList);
 587     return;
 588   }
 589
 590   // run fast jet
 591   // employ setters for these...
 592
 593   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 594   fastjet::AreaType areaType =   fastjet::active_area;
 595   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 596   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 597   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
 598
 599   //range where to compute background
 600   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 601   phiMin = 0;
 602   phiMax = 2*TMath::Pi();
 603   rapMax = fGhostEtamax - fRparam;
 604   rapMin = - fGhostEtamax + fRparam;
 605   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 606
 607
 608   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
 609   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
 610
 611
 612  // loop over all jets and fill information, first one is the leading jet
 613
 614   if(inclusiveJets.size()>0){
 615
 616     //background estimates:all bckg jets wo the 2 hardest
 617     vector<fastjet::PseudoJet> jets2=sortedJets;
 618     if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); //removes the two jets with the highest pT; +2 is correct ro remove 2 jets
 619     Double_t bkg1=0;
 620     Double_t sigma1=0.;
 621     Double_t meanarea1=0.;
 622     clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
 623     background = bkg1;//sets background variable of the task to the correct value
 624
 625
 626     // generate random cones
 627     if(fTCARandomConesOut){
 628       // create a random jet within the acceptance
 629       Double_t etaMax = fTrackEtaWindow - fRparam;//0.9 - 0.4
 630       Int_t nCone = 0;
 631       Int_t nConeRan = 0;
 632       Double_t pTC = 1; // small number
 633       Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
 634       Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
 635       // use fixed position for random Cones
 636       if(randCone_pos){
 637         etaC = randCone_Eta;
 638         phiC = randCone_Phi;
 639       }
 640       // massless jet
 641       Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 642       Double_t pZC = pTC/TMath::Tan(thetaC);
 643       Double_t pXC = pTC * TMath::Cos(phiC);
 644       Double_t pYC = pTC * TMath::Sin(phiC);
 645       Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 646       AliAODJet tmpRecC (pXC,pYC,pZC, pC);
 647
 648       tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
 649       if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
 650       if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
 651
 652       // loop over the reconstructed particles and add up the pT in the random cones
 653       // maybe better to loop over randomized particles not in the real jets...
 654       // but this by definition brings dow average energy in the whole  event
 655       AliAODJet vTmpRanR(1,0,0,1);
 656       for(int i = 0; i < recParticles.GetEntries(); i++){
 657         AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 658         //add up energy in cone
 659         if(fTCARandomConesOut){
 660           AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(0);
 661           if(jC&&jC->DeltaR(vp)<fRparam){
 662             if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
 663             jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
 664           }
 665         }// add up energy in cone
 666
 667         // the randomized input changes eta and phi, but keeps the p_T
 668         if(fTCARandomConesOutRan){
 669           Double_t pTR = vp->Pt();
 670           Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
 671           Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
 672
 673           Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
 674           Double_t pZR = pTR/TMath::Tan(thetaR);
 675
 676           Double_t pXR = pTR * TMath::Cos(phiR);
 677           Double_t pYR = pTR * TMath::Sin(phiR);
 678           Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
 679           vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
 680
 681           AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(0);
 682           if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
 683             if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
 684             jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
 685           }
 686         }
 687       }// loop over recparticles
 688     }  //fTCARandomConesOut
 689     Float_t jetArea = fRparam*fRparam*TMath::Pi();
 690     if(fTCARandomConesOut){
 691       // rescale the momentum vectors for the random cones
 692       AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(0);
 693       if(rC){
 694         Double_t etaC = rC->Eta();
 695         Double_t phiC = rC->Phi();
 696         // massless jet, unit vector
 697         Double_t pTC = rC->ChargedBgEnergy();
 698         if(pTC<=0)pTC = 0.001; // for almost empty events
 699         Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 700         Double_t pZC = pTC/TMath::Tan(thetaC);
 701         Double_t pXC = pTC * TMath::Cos(phiC);
 702         Double_t pYC = pTC * TMath::Sin(phiC);
 703         Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 704         rC->SetPxPyPzE(pXC,pYC,pZC, pC);
 705         rC->SetBgEnergy(0,0);
 706         rC->SetEffArea(jetArea,0);
 707       }
 708     }
 709     if(fTCARandomConesOutRan){
 710       AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(0);
 711       // same with random
 712       if(rC){
 713         Double_t etaC = rC->Eta();
 714         Double_t phiC = rC->Phi();
 715         // massless jet, unit vector
 716         Double_t pTC = rC->ChargedBgEnergy();
 717         if(pTC<=0)pTC = 0.001;// for almost empty events
 718         Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 719         Double_t pZC = pTC/TMath::Tan(thetaC);
 720         Double_t pXC = pTC * TMath::Cos(phiC);
 721         Double_t pYC = pTC * TMath::Sin(phiC);
 722         Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 723         rC->SetPxPyPzE(pXC,pYC,pZC, pC);
 724         rC->SetBgEnergy(0,0);
 725         rC->SetEffArea(jetArea,0);
 726       }
 727     }//find the random jets
 728   }//inclusive Jets > 0
 729
 730  //Calculate delta pT
 731  AliAODJet *randCone = (AliAODJet*)fTCARandomConesOut->At(0);
 732  if(randCone){
 733    //background is the backbround density per area and area=pi*0.4^2 -> backgroundCone is the background energy under the cone
 734    Float_t backgroundCone = background * randCone->EffectiveAreaCharged();
 735    //calculates difference between expected and measured energy density
 736    Float_t ptSub = randCone->Pt() - backgroundCone;
 737    fh1DeltapT->Fill(ptSub);// delta pT
 738    fh1Rho->Fill(background);// background rho
 739    fh1PtRandCone->Fill(randCone->Pt());// pT of random cone
 740    fh1Area->Fill(randCone->EffectiveAreaCharged()); // area of random cone; should always be pi*0.4^2 = 0.5
 741  }else{
 742    if(fDebug)Printf("%s:%d No random Cone found",(char*)__FILE__,__LINE__);
 743  }
 744
 745
 746  if (fDebug > 2){
 747    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
 748    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
 749  }
 750  PostData(1, fHistList);
 751 }
 752
 753 void AliAnalysisTaskJetHBOM::Terminate(Option_t */*option*/)
 754 {
 755   //
 756   // Terminate analysis
 757   //
 758     if (fDebug > 1) printf("AnalysisJetHBOM: Terminate() \n");
 759
 760     if(fMomResH1Fit) delete fMomResH1Fit;
 761     if(fMomResH2Fit) delete fMomResH2Fit;
 762     if(fMomResH3Fit) delete fMomResH3Fit;
 763
 764 }
 765
 766
 767 Int_t  AliAnalysisTaskJetHBOM::GetListOfTracks(TList *list,Int_t type){
 768
 769   //
 770   // get list of tracks/particles for different types
 771   //
 772
 773   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 774
 775   Int_t iCount = 0;
 776    if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
 777     if(type!=kTrackAODextraonly) {
 778       AliAODEvent *aod = 0;
 779       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 780       else aod = AODEvent();
 781       if(!aod){
 782         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 783         return iCount;
 784       }
 785
 786       for(int it = 0;it < aod->GetNumberOfTracks();++it){
 787         AliAODTrack *tr = aod->GetTrack(it);
 788         Bool_t bGood = false;
 789         if(fFilterType == 0)bGood = true;
 790         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
 791         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
 792         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
 793           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
 794           continue;
 795         }
 796         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
 797           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
 798           continue;
 799         }
 800         if(tr->Pt()<fTrackPtCut){
 801           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
 802           continue;
 803         }
 804         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
 805         list->Add(tr);
 806         iCount++;
 807       }
 808     }
 809     if(type==kTrackAODextra || type==kTrackAODextraonly) {
 810       AliAODEvent *aod = 0;
 811       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 812       else aod = AODEvent();
 813
 814       if(!aod){
 815         return iCount;
 816       }
 817       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
 818       if(!aodExtraTracks)return iCount;
 819       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
 820         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
 821         if (!track) continue;
 822
 823         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
 824         if(!trackAOD)continue;
 825         Bool_t bGood = false;
 826         if(fFilterType == 0)bGood = true;
 827         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
 828         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
 829         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
 830         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
 831         if(trackAOD->Pt()<fTrackPtCut) continue;
 832         list->Add(trackAOD);
 833         iCount++;
 834       }
 835     }
 836   }
 837   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
 838     AliMCEvent* mcEvent = MCEvent();
 839     if(!mcEvent)return iCount;
 840     // we want to have alivpartilces so use get track
 841     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
 842       if(!mcEvent->IsPhysicalPrimary(it))continue;
 843       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
 844       if(type == kTrackKineAll){
 845         if(part->Pt()<fTrackPtCut)continue;
 846         list->Add(part);
 847         iCount++;
 848       }
 849       else if(type == kTrackKineCharged){
 850         if(part->Particle()->GetPDG()->Charge()==0)continue;
 851         if(part->Pt()<fTrackPtCut)continue;
 852         list->Add(part);
 853         iCount++;
 854       }
 855     }
 856   }
 857   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 858     AliAODEvent *aod = 0;
 859     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
 860     else aod = AODEvent();
 861     if(!aod)return iCount;
 862     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 863     if(!tca)return iCount;
 864     for(int it = 0;it < tca->GetEntriesFast();++it){
 865       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
 866       if(!part->IsPhysicalPrimary())continue;
 867       if(type == kTrackAODMCAll){
 868         if(part->Pt()<fTrackPtCut)continue;
 869         list->Add(part);
 870         iCount++;
 871       }
 872       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 873         if(part->Charge()==0)continue;
 874         if(part->Pt()<fTrackPtCut)continue;
 875         if(kTrackAODMCCharged){
 876           list->Add(part);
 877         }
 878         else {
 879           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
 880           list->Add(part);
 881         }
 882         iCount++;
 883       }
 884     }
 885   }// AODMCparticle
 886   list->Sort();
 887   return iCount;
 888 }
 889
 890 void AliAnalysisTaskJetHBOM::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
 891
 892   //
 893   // set mom res profiles
 894   //
 895
 896   fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
 897   fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
 898   fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
 899 }
 900
 901 void AliAnalysisTaskJetHBOM:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
 902   //
 903   // set tracking efficiency histos
 904   //
 905
 906   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
 907   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
 908   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
 909 }
 910
 911 Double_t AliAnalysisTaskJetHBOM::GetMomentumSmearing(Int_t cat, Double_t pt) {
 912
 913   //
 914   // Get smearing on generated momentum
 915   //
 916
 917   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
 918
 919   TProfile *fMomRes = 0x0;
 920   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
 921   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
 922   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
 923
 924   if(!fMomRes) {
 925     return 0.;
 926   }
 927
 928
 929   Double_t smear = 0.;
 930
 931   if(pt>20.) {
 932     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
 933     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
 934     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
 935   }
 936   else {
 937
 938     Int_t bin = fMomRes->FindBin(pt);
 939
 940     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
 941
 942   }
 943
 944   if(fMomRes) delete fMomRes;
 945
 946   return smear;
 947 }
 948
 949 void AliAnalysisTaskJetHBOM::FitMomentumResolution() {
 950   //
 951   // Fit linear function on momentum resolution at high pT
 952   //
 953
 954   if(!fMomResH1Fit && fMomResH1) {
 955     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
 956     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
 957     fMomResH1Fit ->SetRange(5.,100.);
 958   }
 959
 960   if(!fMomResH2Fit && fMomResH2) {
 961     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
 962     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
 963     fMomResH2Fit ->SetRange(5.,100.);
 964   }
 965
 966   if(!fMomResH3Fit && fMomResH3) {
 967     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
 968     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
 969     fMomResH3Fit ->SetRange(5.,100.);
 970   }
 971
 972 }
 973