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