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