JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx

   1 // **************************************
   2 // Task 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 <TH1F.h>
  25 #include <TH2F.h>
  26 #include <THnSparse.h>
  27 #include <TSystem.h>
  28 #include <TInterpreter.h>
  29 #include <TList.h>
  30 #include <TLorentzVector.h>
  31 #include <TRefArray.h>
  32 #include  "TDatabasePDG.h"
  33
  34 #include "AliAnalysisTaskJetBackgroundSubtract.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAODTrack.h"
  38 #include "AliAODJet.h"
  39 #include "AliAODEvent.h"
  40 #include "AliInputEventHandler.h"
  41 #include "AliAODJetEventBackground.h"
  42
  43
  44 ClassImp(AliAnalysisTaskJetBackgroundSubtract)
  45
  46 AliAnalysisTaskJetBackgroundSubtract::~AliAnalysisTaskJetBackgroundSubtract(){
  47   delete fJBArray;
  48   delete fOutJetArrayList;
  49   delete fInJetArrayList;
  50 }
  51
  52 AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract():
  53   AliAnalysisTaskSE(),
  54   fAODOut(0x0),
  55   fAODIn(0x0),
  56   fAODExtension(0x0),
  57   fJBArray(new TObjArray()),
  58   fBackgroundBranch(""),
  59   fNonStdFile(""),
  60   fReplaceString1("B0"),
  61   fReplaceString2("B%d"),
  62   fSubtraction(k4Area),
  63   fKeepJets(kFALSE),
  64   fInJetArrayList(0x0),
  65   fOutJetArrayList(0x0),
  66   fh2CentvsRho(0x0),
  67   fh2CentvsSigma(0x0),
  68   fh2MultvsRho(0x0),
  69   fh2MultvsSigma(0x0),
  70   fh2ShiftEta(0x0),
  71   fh2ShiftPhi(0x0),
  72   fh2ShiftEtaLeading(0x0),
  73   fh2ShiftPhiLeading(0x0),
  74   fHistList(0x0)
  75 {
  76
  77 }
  78
  79 AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract(const char* name):
  80
  81   AliAnalysisTaskSE(name),
  82   fAODOut(0x0),
  83   fAODIn(0x0),
  84   fAODExtension(0x0),
  85   fJBArray(new TObjArray()),
  86   fBackgroundBranch(""),
  87   fNonStdFile(""),
  88   fReplaceString1("B0"),
  89   fReplaceString2("B%d"),
  90   fSubtraction(k4Area),
  91   fKeepJets(kFALSE),
  92   fInJetArrayList(0x0),
  93   fOutJetArrayList(0x0),
  94   fh2CentvsRho(0x0),
  95   fh2CentvsSigma(0x0),
  96   fh2MultvsRho(0x0),
  97   fh2MultvsSigma(0x0),
  98   fh2ShiftEta(0x0),
  99   fh2ShiftPhi(0x0),
 100   fh2ShiftEtaLeading(0x0),
 101   fh2ShiftPhiLeading(0x0),
 102   fHistList(0x0)
 103 {
 104  DefineOutput(1, TList::Class());
 105 }
 106
 107
 108
 109 Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify()
 110 {
 111   //
 112   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 113
 114   ResetOutJets();
 115
 116   // Now we also have the Input Event Available! Fillvthe pointers in the list
 117
 118   fInJetArrayList->Clear();
 119   fOutJetArrayList->Clear();
 120
 121   for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
 122     TObjString *ostr = (TObjString*)fJBArray->At(iJB);
 123
 124
 125     TClonesArray* jarray = 0;
 126     if(!jarray&&fAODOut){
 127       jarray = (TClonesArray*)(fAODOut->FindListObject(ostr->GetString().Data()));
 128     }
 129     if(!jarray&&fAODExtension){
 130       jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(ostr->GetString().Data()));
 131     }
 132     if(!jarray&&fAODIn){
 133       jarray = (TClonesArray*)(fAODIn->FindListObject(ostr->GetString().Data()));
 134     }
 135
 136     if(!jarray){
 137       if(fDebug){
 138         Printf("%s:%d Input jet branch %s not found",(char*)__FILE__,__LINE__,ostr->GetString().Data());
 139       }
 140       continue;
 141     }
 142
 143     TString newName(ostr->GetString().Data());
 144     newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
 145     TClonesArray* jarrayOut = 0;
 146     if(newName.CompareTo(ostr->GetString())==0){
 147       Printf("%s:%d Input and output branch would have the same name, skipping %s ",(char*)__FILE__,__LINE__,ostr->GetString().Data());
 148       continue;
 149     }
 150
 151     if(!jarrayOut&&fAODOut){
 152       jarrayOut = (TClonesArray*)(fAODOut->FindListObject(newName.Data()));
 153     }
 154     if(!jarrayOut&&fAODExtension){
 155       jarrayOut = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(newName.Data()));
 156     }
 157
 158     if(!jarrayOut){
 159       if(fDebug){
 160         Printf("%s:%d Output jet branch %s not found",(char*)__FILE__,__LINE__,newName.Data());
 161         PrintAODContents();
 162       }
 163       continue;
 164     }
 165     if(jarrayOut&&jarray){
 166       fOutJetArrayList->Add(jarrayOut);
 167       fInJetArrayList->Add(jarray);
 168     }
 169   }
 170   return kTRUE;
 171 }
 172
 173 void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
 174 {
 175
 176   //
 177   // Create the output container
 178   //
 179   // Connect the AOD
 180
 181   if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n");
 182
 183
 184
 185   if(fNonStdFile.Length()!=0){
 186
 187     // case that we have an AOD extension we need to fetch the jets from the extended output
 188     // we identifay the extension aod event by looking for the branchname
 189     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 190
 191     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 192
 193     if(!fAODExtension){
 194       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 195     }
 196   }
 197   fAODOut = AODEvent();
 198
 199   // usually we do not have the input already here
 200
 201   if(!fInJetArrayList)fInJetArrayList =new TList();
 202   if(!fOutJetArrayList)fOutJetArrayList =new TList();
 203
 204   for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
 205     TObjString *ostr = (TObjString*)fJBArray->At(iJB);
 206     TString newName(ostr->GetString().Data());
 207     if(!newName.Contains(fReplaceString1.Data())){
 208       Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
 209              newName.Data());
 210       continue;
 211     }
 212
 213
 214     // add a new branch to the output for the background subtracted jets take the names from
 215     // the input jets and replace the background flag names
 216     TClonesArray *tca = new TClonesArray("AliAODJet", 0);
 217     newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
 218     if(newName.CompareTo(ostr->GetString())==0){
 219       Printf("%s:%d Input and output branch would have the same name, skipping: %s ",(char*)__FILE__,__LINE__,ostr->GetString().Data());
 220       continue;
 221     }
 222
 223     if(fDebug){
 224       Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(),
 225              ostr->GetString().Data());
 226     }
 227     tca->SetName(newName.Data());
 228     AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
 229   }
 230
 231
 232   if(!fHistList)fHistList = new TList();
 233   fHistList->SetOwner();
 234   PostData(1, fHistList); // post data in any case once
 235
 236   //
 237
 238   // delta pT vs. area vs. cent vs. mult
 239   const Int_t nSparseBinsDelta = 4;
 240   const Int_t nBinsDelta[nSparseBinsDelta] =   {   241,  10,  10,  25};
 241   const Double_t xminDelta[nSparseBinsDelta] = {-120.5,   0,   0,   0};
 242   const Double_t xmaxDelta[nSparseBinsDelta] = { 120.5, 1.0, 100,5000};
 243
 244   for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
 245     TObjString *ostr = (TObjString*)fJBArray->At(iJB);
 246     TString oldName(ostr->GetString().Data());
 247     TString newName(ostr->GetString().Data());
 248     if(!newName.Contains(fReplaceString1.Data())){
 249       Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
 250              newName.Data());
 251       continue;
 252     }
 253     newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
 254
 255     TH2F *hTmp = new TH2F(Form("h2PtInPtOut_%d",iJB),Form(";%s p_{T}; %s p_{T}",oldName.Data(),newName.Data()),200,0,200.,400,-200.,200.);
 256     fHistList->Add(hTmp);
 257     THnSparseF *hFTmp = new THnSparseF(Form("hnDPtAreaCentMult_%d",iJB),Form("%s delta;#delta p_{T};Area;cent;mult",newName.Data()),nSparseBinsDelta,nBinsDelta,xminDelta,xmaxDelta);
 258     fHistList->Add(hFTmp);
 259   }
 260
 261   Bool_t oldStatus = TH1::AddDirectoryStatus();
 262   TH1::AddDirectory(kFALSE);
 263
 264   //
 265   //  Histogram booking, add som control histograms here
 266   //
 267
 268
 269     fh2CentvsRho = new TH2F("fh2CentvsRho","centrality vs background density", 100,0.,100.,600,0.,300.);
 270     fh2CentvsSigma = new TH2F("fh2CentvsSigma","centrality vs backgroun sigma",100,0.,100.,500,0.,50.);
 271     fHistList->Add(fh2CentvsRho);
 272     fHistList->Add(fh2CentvsSigma);
 273
 274     fh2MultvsRho = new TH2F("fh2MultvsRho","mult vs background density", 100,0.,5000.,600,0.,300.);
 275     fh2MultvsSigma = new TH2F("fh2MultvsSigma","mult vs background sigma",100,0.,5000.,500,0.,50.);
 276     fHistList->Add(fh2MultvsRho);
 277     fHistList->Add(fh2MultvsSigma);
 278
 279    if(fSubtraction==k4Area){
 280     fh2ShiftEta = new TH2F("fh2ShiftEta","extended correction Eta",100,-0.9,0.9,100,-0.9,0.9);
 281     fh2ShiftPhi = new TH2F("fh2ShiftPhi","extended correction Phi",100,0.,6.5,100,0.,6.5);
 282     fh2ShiftEtaLeading = new TH2F("fh2ShiftEtaLeading","extended correction Eta",100,-0.9,0.9,100,-0.9,0.9);
 283     fh2ShiftPhiLeading = new TH2F("fh2ShiftPhiLeading","extended correction Phi",100,0.,6.5,100,0.,6.5);
 284      fHistList->Add(fh2ShiftEta);
 285      fHistList->Add(fh2ShiftPhi);
 286      fHistList->Add(fh2ShiftEtaLeading);
 287      fHistList->Add(fh2ShiftPhiLeading);
 288    }
 289
 290   // =========== Switch on Sumw2 for all histos ===========
 291   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 292     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 293     if (h1){
 294       h1->Sumw2();
 295       continue;
 296     }
 297     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 298     if(hn)hn->Sumw2();
 299   }
 300   TH1::AddDirectory(oldStatus);
 301
 302   if(fBackgroundBranch.Length()==0)
 303     AliError(Form("%s:%d No BackgroundBranch defined",(char*)__FILE__,__LINE__));
 304   if(fJBArray->GetEntries()==0)
 305     AliError(Form("%s:%d No Jet Branches defined defined",(char*)__FILE__,__LINE__));
 306 }
 307
 308 void AliAnalysisTaskJetBackgroundSubtract::Init()
 309 {
 310   //
 311   // Initialization
 312   //
 313   if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::Init() \n");
 314 }
 315
 316 void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
 317 {
 318
 319   if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
 320   ResetOutJets();
 321   if(fBackgroundBranch.Length()==0||fJBArray->GetEntries()==0){
 322     if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
 323     PostData(1,fHistList);
 324   }
 325   if(fJBArray->GetEntries()!=fInJetArrayList->GetEntries()){
 326     if(fDebug)Printf("%s:%d different Array  sizes %d %d %d",(char*)__FILE__,__LINE__,fJBArray->GetEntries(),fInJetArrayList->GetEntries(),fOutJetArrayList->GetEntries());
 327     PostData(1,fHistList);
 328   }
 329
 330
 331
 332   AliAODJetEventBackground*  evBkg = 0;
 333   TClonesArray*              bkgClusters = 0;
 334   TClonesArray*              bkgClustersRC = 0;
 335   TString bkgClusterName(fBackgroundBranch.Data());
 336   bkgClusterName.ReplaceAll(Form("%s_",AliAODJetEventBackground::StdBranchName()),"");
 337   TString bkgClusterRCName(Form("%s%s",bkgClusterName.Data(),"RandomCone"));
 338
 339   if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODOut){
 340     evBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
 341     bkgClusters = (TClonesArray*)(fAODOut->FindListObject(bkgClusterName.Data()));
 342     bkgClustersRC = (TClonesArray*)(fAODOut->FindListObject(bkgClusterRCName.Data()));
 343
 344     if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
 345     if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
 346     if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
 347   }
 348   if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODExtension){
 349     evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 350     bkgClusters = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterName.Data()));
 351     bkgClustersRC = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterRCName.Data()));
 352     if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
 353     if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
 354
 355     if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
 356   }
 357
 358   if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODIn){
 359     evBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
 360     bkgClusters = (TClonesArray*)(fAODIn->FindListObject(bkgClusterName.Data()));
 361     bkgClustersRC = (TClonesArray*)(fAODIn->FindListObject(bkgClusterRCName.Data()));
 362
 363     if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
 364     if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
 365     if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
 366   }
 367
 368   if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc||fSubtraction==k4Area)){
 369     if(fDebug){
 370       Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
 371       PrintAODContents();
 372     }
 373     PostData(1,fHistList);
 374     return;
 375   }
 376
 377   if(!bkgClusters&&(fSubtraction==kRhoRecalc)){
 378     if(fDebug){
 379       Printf("%s:%d Background cluster branch %s not found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
 380       PrintAODContents();
 381     }
 382     PostData(1,fHistList);
 383     return;
 384   }
 385
 386   if(!bkgClustersRC&&(fSubtraction==kRhoRC)){
 387     if(fDebug){
 388       Printf("%s:%d Background cluster RC branch %s not found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
 389       PrintAODContents();
 390     }
 391     PostData(1,fHistList);
 392     return;
 393   }
 394   // LOOP over all jet branches and subtract the background
 395
 396    Float_t rho = 0;
 397    Float_t sigma=0.;
 398    Double_t meanarea = 0;
 399    TLorentzVector backgroundv;
 400    Float_t cent=0.;
 401
 402    if(fAODOut)cent = fAODOut->GetHeader()->GetCentrality();
 403    if(fAODIn) cent = fAODIn->GetHeader()->GetCentrality();
 404
 405    if(evBkg)sigma=evBkg->GetSigma(1);
 406
 407    if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
 408    if(fSubtraction==k4Area){
 409      rho = evBkg->GetBackground(0);
 410      sigma=evBkg->GetSigma(0);
 411    }
 412
 413    if(fSubtraction==kRhoRecalc){
 414      meanarea=evBkg->GetMeanarea(1);
 415      rho =RecalcRho(bkgClusters,meanarea);
 416    }
 417    if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
 418
 419    Float_t mult = 0;
 420    for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
 421     TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
 422     if(jarray){
 423       TString tmp(jarray->GetName());
 424       if(tmp.Contains("cluster")){
 425         mult = MultFromJetRefs(jarray);
 426         if(mult>0)break;
 427       }
 428     }
 429    }
 430
 431    fh2CentvsRho->Fill(cent,rho);
 432    fh2CentvsSigma->Fill(cent,sigma);
 433
 434    fh2MultvsRho->Fill(mult,rho);
 435    fh2MultvsSigma->Fill(mult,sigma);
 436
 437    for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
 438     TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
 439     TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
 440
 441     if(!jarray||!jarrayOut){
 442       Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray,jarrayOut);
 443       continue;
 444     }
 445     TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
 446     THnSparseF* hnDPtAreaCentMult = (THnSparseF*)fHistList->FindObject(Form("hnDPtAreaCentMult_%d",iJB));
 447     // loop over all jets
 448     Int_t nOut = 0;
 449
 450     Double_t deltaPt[4];
 451     deltaPt[2] = cent;
 452     deltaPt[3] = mult;
 453
 454     for(int i = 0;i < jarray->GetEntriesFast();i++){
 455       AliAODJet *jet = (AliAODJet*)jarray->At(i);
 456       AliAODJet tmpNewJet(*jet);
 457       Bool_t bAdd = false;
 458       Float_t ptSub = 0;
 459
 460
 461       if(fSubtraction==kArea){
 462         Double_t background = rho * jet->EffectiveAreaCharged();
 463         ptSub = jet->Pt() - background;
 464         if(fDebug>2){
 465           Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
 466         }
 467         if(ptSub<=0){
 468           // optionally rescale it and keep??
 469           if(fKeepJets){
 470              bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
 471           }
 472           else{
 473             bAdd = false;
 474           }
 475         }
 476         else{
 477           bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
 478         }
 479         // add background estimates to the new jet object
 480         // allows to recover old p_T and rho...
 481         tmpNewJet.SetBgEnergy(background,0);
 482         tmpNewJet.SetPtSubtracted(ptSub,0);
 483       }// kAREA
 484       else if(fSubtraction==kRhoRecalc){
 485         Double_t background = rho * jet->EffectiveAreaCharged();
 486         ptSub = jet->Pt() - background;
 487         if(fDebug>2){
 488           Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
 489         if(ptSub<=0){
 490           // optionally rescale it and keep
 491           if(fKeepJets){
 492              bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
 493           }
 494           else{
 495             bAdd = false;
 496           }
 497         }
 498         else{
 499           bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
 500         }
 501         // add background estimates to the new jet object
 502         // allows to recover old p_T and rho...
 503         tmpNewJet.SetBgEnergy(background,0);
 504         tmpNewJet.SetPtSubtracted(ptSub,0);
 505       }//kRhoRecalc
 506        else if(fSubtraction==kRhoRC){
 507         Double_t background = rho * jet->EffectiveAreaCharged();
 508         ptSub = jet->Pt() - background;
 509         if(fDebug>2){   Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
 510         if(ptSub<=0){
 511           if(fKeepJets){
 512              bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
 513           }
 514           else{
 515             bAdd = false;
 516           }
 517         }
 518         else{
 519           bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
 520         }
 521         // add background estimates to the new jet object
 522         // allows to recover old p_T and rho...
 523         tmpNewJet.SetBgEnergy(background,0);
 524         tmpNewJet.SetPtSubtracted(ptSub,0);
 525        }//kRhoRC
 526
 527        else if(fSubtraction==k4Area&&jet->VectorAreaCharged()){
 528          backgroundv.SetPxPyPzE(rho*(jet->VectorAreaCharged())->Px(),rho*(jet->VectorAreaCharged())->Py(),rho*(jet->VectorAreaCharged())->Pz(),rho*(jet->VectorAreaCharged())->E());
 529          ptSub = jet->Pt()-backgroundv.Pt();
 530          if((backgroundv.E()>=jet->E())||(backgroundv.Pt()>=jet->Pt())){
 531            if(fKeepJets){
 532              bAdd =  RescaleJetMomentum(&tmpNewJet,0.1);
 533            }
 534            else{
 535              bAdd = false;
 536            }
 537          }
 538          else{
 539            bAdd = RescaleJet4vector(&tmpNewJet,backgroundv);
 540          }
 541          // add background estimates to the new jet object
 542          // allows to recover old p_T and rho...
 543          tmpNewJet.SetBgEnergy(backgroundv.Pt(),0);
 544          tmpNewJet.SetPtSubtracted(ptSub,0);
 545
 546        }//kArea4vector
 547
 548       if(bAdd){
 549         AliAODJet *newJet = new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
 550         // what about track references, clear for now...
 551         if(fSubtraction==k4Area){
 552          fh2ShiftEta->Fill(jet->Eta(),newJet->Eta());
 553          fh2ShiftPhi->Fill(jet->Phi(),newJet->Phi());
 554          if(i==0){fh2ShiftEtaLeading->Fill(jet->Eta(),newJet->Eta());
 555            fh2ShiftPhiLeading->Fill(jet->Phi(),newJet->Phi());}}
 556
 557         // set the references
 558         newJet->GetRefTracks()->Clear();
 559         TRefArray *refs = jet->GetRefTracks();
 560         for(Int_t ir=0;ir<refs->GetEntriesFast();ir++){
 561           AliVParticle *vp = dynamic_cast<AliVParticle*>(refs->At(ir));
 562           if(vp)newJet->AddTrack(vp);
 563         }
 564       }
 565       if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
 566       if(hnDPtAreaCentMult){
 567         deltaPt[0] = ptSub;
 568         deltaPt[1] = jet->EffectiveAreaCharged();
 569         hnDPtAreaCentMult->Fill(deltaPt);
 570       }
 571     }
 572     if(jarrayOut)jarrayOut->Sort();
 573    }
 574
 575    PostData(1, fHistList);
 576 }
 577
 578 void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t */*option*/)
 579 {
 580   // Terminate analysis
 581   //
 582   if (fDebug > 1) printf("AnalysisJetBackgroundSubtract: Terminate() \n");
 583 }
 584
 585 Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJetMomentum(AliAODJet *jet,Float_t pT){
 586   // keep the direction and the jet mass
 587   if(pT<=0)return kFALSE;
 588   Double_t pTold = jet->Pt();
 589   Double_t scale  = pT/pTold;
 590   Double_t mass  = jet->M();
 591   Double_t pNew = jet->P() * scale;
 592   jet->SetPxPyPzE(scale*jet->Px(),scale*jet->Py(),scale*jet->Pz(),TMath::Sqrt(mass*mass+pNew*pNew));
 593
 594
 595
 596   return kTRUE;
 597 }
 598
 599 Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJet4vector(AliAODJet *jet,TLorentzVector backgroundv){
 600
 601   if(backgroundv.Pt()<0.) return kFALSE;
 602   jet->SetPxPyPzE(jet->Px()-backgroundv.Px(),jet->Py()-backgroundv.Py(),jet->Pz()-backgroundv.Pz(),jet->E()-backgroundv.E());
 603
 604  return kTRUE;
 605 }
 606
 607
 608
 609
 610
 611
 612
 613
 614 Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
 615
 616   Double_t ptarea=0.;
 617   Int_t count=0;
 618   Double_t rho=0.;
 619   const Double_t Rlimit2=0.8*0.8;  //2*jet radius.
 620   TClonesArray* jarray=0;
 621
 622   for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
 623     TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
 624     TString jetref=ostr->GetString().Data();
 625     if(jetref.Contains("ANTIKT04")){
 626       jarray = (TClonesArray*)fInJetArrayList->At(iJB);
 627     }
 628   }
 629   if(!jarray)return rho;
 630   if(jarray->GetEntries()>=2){
 631     AliAODJet *first = (AliAODJet*)(jarray->At(0));
 632     AliAODJet *second= (AliAODJet*)(jarray->At(1));
 633     for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
 634       AliAODJet *clus = (AliAODJet*)(bkgClusters->At(k));
 635       if(TMath::Abs(clus->Eta())>0.5) continue;
 636       if((clus->EffectiveAreaCharged())<0.1*meanarea) continue;
 637       Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
 638         (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
 639       Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
 640         (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
 641       if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;
 642       ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
 643       count=count+1;}
 644     if(count!=0) rho=ptarea/count;
 645   }
 646   return rho;
 647 }
 648
 649    Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
 650
 651        Double_t ptarea=0.;
 652        Int_t count=0;
 653        Double_t rho=0.;
 654        const Double_t Rlimit2=0.8*0.8;  //2*jet radius.
 655        TClonesArray* jarray=0;
 656         for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
 657           TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
 658           TString jetref=ostr->GetString().Data();
 659           if(jetref.Contains("ANTIKT04")){
 660             jarray = (TClonesArray*)fInJetArrayList->At(iJB);
 661           }
 662         }
 663         if(!jarray)return rho;
 664
 665          if(jarray->GetEntries()>=2){
 666            AliAODJet *first = (AliAODJet*)(jarray->At(0));
 667            AliAODJet *second=(AliAODJet*)(jarray->At(1));
 668          for(Int_t k=0;k<bkgClustersRC->GetEntriesFast();k++){
 669            AliAODJet *clus = (AliAODJet*)(bkgClustersRC->At(k));
 670            if(TMath::Abs(clus->Eta())>0.5) continue;
 671            Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
 672              (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
 673            Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
 674              (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
 675            if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;
 676            ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
 677            count=count+1;}
 678          if(count!=0) rho=ptarea/count;  }
 679          return rho;
 680 }
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690 void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
 691   if(!fOutJetArrayList)return;
 692   for(int iJB = 0;iJB<fOutJetArrayList->GetEntries();iJB++){
 693     TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB);
 694     if(jarray)jarray->Delete();
 695   }
 696 }
 697
 698
 699 void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
 700   if(fAODIn){
 701     Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
 702     fAODIn->Print();
 703   }
 704   if(fAODExtension){
 705     Printf("%s:%d >>>>>> Extenstion",(char*)__FILE__,__LINE__);
 706     fAODExtension->GetAOD()->Print();
 707   }
 708   if(fAODOut){
 709     Printf("%s:%d >>>>>> Output",(char*)__FILE__,__LINE__);
 710     fAODOut->Print();
 711   }
 712 }
 713
 714 Int_t AliAnalysisTaskJetBackgroundSubtract::MultFromJetRefs(TClonesArray* jets){
 715   if(!jets)return 0;
 716
 717   Int_t refMult = 0;
 718   for(int ij = 0;ij < jets->GetEntries();++ij){
 719     AliAODJet* jet = (AliAODJet*)jets->At(ij);
 720     if(!jet)continue;
 721     TRefArray *refs = jet->GetRefTracks();
 722     if(!refs)continue;
 723     refMult += refs->GetEntries();
 724   }
 725   return refMult;
 726
 727 }