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