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