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