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