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