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