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