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