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