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