Adding includes for EMCAL_Utils and OADB PATH (A. Shabetai)
[u/mrichter/AliRoot.git] / 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         // set the references 
497         newJet->GetRefTracks()->Clear();
498         TRefArray *refs = jet->GetRefTracks();
499         for(Int_t ir=0;ir<refs->GetEntriesFast();ir++){
500           AliVParticle *vp = dynamic_cast<AliVParticle*>(refs->At(ir));
501           if(vp)newJet->AddTrack(vp);
502         }
503       }
504     }
505     if(jarrayOut)jarrayOut->Sort();
506    }
507    
508    PostData(1, fHistList);
509 }
510
511 void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t */*option*/)
512 {
513   // Terminate analysis
514   //
515   if (fDebug > 1) printf("AnalysisJetBackgroundSubtract: Terminate() \n");
516 }
517
518 Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJetMomentum(AliAODJet *jet,Float_t pT){
519   // keep the direction and the jet mass
520   if(pT<=0)return kFALSE;
521   Double_t pTold = jet->Pt();
522   Double_t scale  = pT/pTold;
523   Double_t mass  = jet->M();
524   Double_t pNew = jet->P() * scale;
525   jet->SetPxPyPzE(scale*jet->Px(),scale*jet->Py(),scale*jet->Pz(),TMath::Sqrt(mass*mass+pNew*pNew));
526  
527
528
529   return kTRUE;
530 }
531
532 Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJet4vector(AliAODJet *jet,TLorentzVector backgroundv){
533   
534   if(backgroundv.Pt()<0.) return kFALSE;
535   jet->SetPxPyPzE(jet->Px()-backgroundv.Px(),jet->Py()-backgroundv.Py(),jet->Pz()-backgroundv.Pz(),jet->E()-backgroundv.E());
536  
537  return kTRUE;
538 }
539
540
541
542
543
544
545
546
547 Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
548   
549   Double_t ptarea=0.;
550   Int_t count=0;
551   Double_t rho=0.; 
552   const Double_t Rlimit2=0.8*0.8;  //2*jet radius.
553   TClonesArray* jarray=0;
554   
555   for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
556     TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
557     TString jetref=ostr->GetString().Data();
558     if(jetref.Contains("ANTIKT04")){ 
559       jarray = (TClonesArray*)fInJetArrayList->At(iJB);
560     }
561   }
562   if(!jarray)return rho;
563   if(jarray->GetEntries()>=2){ 
564     AliAODJet *first = (AliAODJet*)(jarray->At(0)); 
565     AliAODJet *second= (AliAODJet*)(jarray->At(1)); 
566     for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
567       AliAODJet *clus = (AliAODJet*)(bkgClusters->At(k));
568       if(TMath::Abs(clus->Eta())>0.5) continue;
569       if((clus->EffectiveAreaCharged())<0.1*meanarea) continue; 
570       Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
571         (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
572       Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
573         (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
574       if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;    
575       ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged(); 
576       count=count+1;}
577     if(count!=0) rho=ptarea/count; 
578   }        
579   return rho;
580 }
581
582    Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
583   
584        Double_t ptarea=0.;
585        Int_t count=0;
586        Double_t rho=0.; 
587        const Double_t Rlimit2=0.8*0.8;  //2*jet radius.
588        TClonesArray* jarray=0;
589         for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
590           TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
591           TString jetref=ostr->GetString().Data();
592           if(jetref.Contains("ANTIKT04")){ 
593             jarray = (TClonesArray*)fInJetArrayList->At(iJB);
594           }
595         }
596         if(!jarray)return rho;
597
598          if(jarray->GetEntries()>=2){ 
599            AliAODJet *first = (AliAODJet*)(jarray->At(0)); 
600            AliAODJet *second=(AliAODJet*)(jarray->At(1)); 
601          for(Int_t k=0;k<bkgClustersRC->GetEntriesFast();k++){
602            AliAODJet *clus = (AliAODJet*)(bkgClustersRC->At(k));
603            if(TMath::Abs(clus->Eta())>0.5) continue;
604            Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
605              (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
606            Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
607              (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
608            if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;    
609            ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged(); 
610            count=count+1;}
611          if(count!=0) rho=ptarea/count;  }
612          return rho;
613 }
614
615
616
617
618
619
620
621
622
623 void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
624   if(!fOutJetArrayList)return;
625   for(int iJB = 0;iJB<fOutJetArrayList->GetEntries();iJB++){
626     TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB);
627     if(jarray)jarray->Delete();
628   }
629 }
630
631
632 void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
633   if(fAODIn){
634     Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
635     fAODIn->Print();
636   }
637   if(fAODExtension){
638     Printf("%s:%d >>>>>> Extenstion",(char*)__FILE__,__LINE__);
639     fAODExtension->GetAOD()->Print();
640   }
641   if(fAODOut){
642     Printf("%s:%d >>>>>> Output",(char*)__FILE__,__LINE__);
643     fAODOut->Print();
644   }
645 }