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