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