Adding the possibility not to follow particle mothers during StepManager (default...
[u/mrichter/AliRoot.git] / JETAN / AliJetBkg.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //--------------------------------------------------
19 // Method implementation for background studies and background subtraction with UA1 algorithms
20 //
21 // Author: magali.estienne@subatech.in2p3.fr
22 //-------------------------------------------------
23
24 #include <Riostream.h> 
25 #include <TList.h>
26 #include <TClonesArray.h>
27 #include <TH1F.h>
28 #include <TH2F.h>
29
30 #include "AliAODJetEventBackground.h"
31 #include "AliUA1JetHeaderV1.h"
32 #include "AliJetCalTrk.h"
33 #include "AliJetBkg.h"
34
35
36 using namespace std;
37
38 ClassImp(AliJetBkg)
39
40 ////////////////////////////////////////////////////////////////////////
41
42 AliJetBkg::AliJetBkg():
43   TObject(),
44   fEvent(0x0),
45   fHeader(0x0),
46   fDebug(0),
47   fhEtBackg(0x0),
48   fhAreaBackg(0x0)
49 {
50   // Default constructor
51   for(int i = 0;i < kMaxJets;i++){
52     fhAreaJet[i] = fhEtJet[i] = 0;
53   }
54 }
55
56 //----------------------------------------------------------------
57 AliJetBkg::AliJetBkg(const AliJetBkg& input):
58   TObject(input),
59   fEvent(input.fEvent),
60   fHeader(input.fHeader),
61   fDebug(input.fDebug),
62   fhEtBackg(input.fhEtBackg),
63   fhAreaBackg(input.fhAreaBackg)
64 {
65   // copy constructor
66   for(int i = 0;i < kMaxJets;i++){
67     fhAreaJet[i] = input.fhAreaJet[i];
68     fhEtJet[i] = input.fhEtJet[i];
69   }
70
71 }
72
73 //----------------------------------------------------------------
74 AliJetBkg::~AliJetBkg()
75 {
76   // Destructor
77   if(fhEtBackg) delete  fhEtBackg;
78   if(fhAreaBackg) delete  fhAreaBackg;
79    for(int i = 0;i < kMaxJets;i++){
80      if(fhAreaJet[i]) delete fhAreaJet[i];
81      if(fhEtJet[i])  delete fhEtJet[i];
82    }
83
84 }
85
86 //----------------------------------------------------------------
87 Bool_t AliJetBkg::PtCutPass(Int_t id, Int_t nTracks)
88 {
89   // Check if track or cell passes the cut flag
90   if(id < nTracks && fEvent->GetCalTrkTrack(id)->GetCutFlag() == 1) 
91    return kTRUE;
92   else return kFALSE;
93
94 }
95
96 //----------------------------------------------------------------
97 Bool_t AliJetBkg::SignalCutPass(Int_t id, Int_t nTracks)
98 {
99   // Check if track or cell passes the cut flag
100   if(id < nTracks && fEvent->GetCalTrkTrack(id)->GetSignalFlag() == 1)
101     return kTRUE;
102   else return kFALSE;
103
104 }
105
106 //----------------------------------------------------------------
107 Float_t AliJetBkg::CalcJetAreaEtaCut(Float_t radius, Float_t etaJet)
108 {
109   // Calculate jet area taking into account an acceptance cut in eta
110   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
111   Float_t detamax = etaJet + radius;
112   Float_t detamin = etaJet - radius;
113   Float_t accmax = 0.0; Float_t accmin = 0.0;
114   if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
115     Float_t h = header->GetLegoEtaMax() - etaJet;
116     accmax = radius*radius*TMath::ACos(h/radius) - h*TMath::Sqrt(radius*radius - h*h);
117   }
118   if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
119     Float_t h = header->GetLegoEtaMax() + etaJet;
120     accmin = radius*radius*TMath::ACos(h/radius) - h*TMath::Sqrt(radius*radius - h*h);
121   }
122   
123   return radius*radius*TMath::Pi() - accmax - accmin;
124
125 }
126
127 //----------------------------------------------------------------
128 void AliJetBkg::CalcJetAndBckgAreaEtaCut(Bool_t calcOutsideArea, Float_t radius, Int_t nJets, const Float_t* etaJet, Float_t* &areaJet, Float_t &areaOut)
129 {
130   // Calculate jet and bacground areas taking into account an acceptance cut in eta
131
132   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
133   areaOut = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax() - header->GetLegoPhiMin());
134   for(Int_t k=0; k<nJets; k++){
135     areaJet[k] = CalcJetAreaEtaCut(radius, etaJet[k]);
136     if(calcOutsideArea) areaOut = areaOut - areaJet[k];
137   }
138
139 }
140
141 //----------------------------------------------------------------
142 void AliJetBkg::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, Float_t&sigmaN,
143                               const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
144                               Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
145                               Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet,
146                               Int_t* injet, Float_t* &areaJet)
147 {
148   //
149   // Background subtraction using cone method but without correction in dE/deta distribution
150   // Cases to take into account the EMCal geometry are included
151   //
152   
153   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
154   //calculate energy inside and outside cones
155   fDebug = header->GetDebug();
156   Float_t rc = header->GetRadius();
157   Float_t etOut = 0;
158   // Get number of tracks from EventCalTrk
159   Int_t nTracks = fEvent->GetNCalTrkTracks();
160
161   Float_t etIn[kMaxJets] = {0};
162   Float_t areaOut = 0.;
163
164   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
165
166     for(Int_t ijet=0; ijet<nJ; ijet++){
167       Float_t deta = etaT[jpart] - etaJet[ijet];
168       Float_t dphi = phiT[jpart] - phiJet[ijet];
169       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
170       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
171       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
172       if(dr <= rc){ // particles inside this cone
173         if(jpart < nTracks) multJetT[ijet]++;
174         else multJetC[ijet]++;
175         multJet[ijet]++;
176         injet[jpart] = ijet;
177         if(PtCutPass(jpart,nTracks)){ // pt cut
178           etIn[ijet] += ptT[jpart];
179           if(SignalCutPass(jpart,nTracks))
180             etsigJet[ijet]+= ptT[jpart];
181         }
182         break;
183       }
184     }// end jets loop
185
186     if((injet[jpart] == -1) &&
187        (PtCutPass(jpart,nTracks))){
188       etOut += ptT[jpart]; // particle outside cones and pt cut
189     }
190   } //end particle loop
191
192   // Calculate jet and background areas
193   Bool_t calcAreaOut = kTRUE;
194   CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut);
195
196   //subtract background using area method
197   for(Int_t ljet=0; ljet<nJ; ljet++){
198     Float_t areaRatio = areaJet[ljet]/areaOut;
199     etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
200   }
201
202   // estimate new total background
203   Float_t areaT = 0;
204   areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
205   etbgTotalN = etOut*areaT/areaOut;
206
207   // estimate standard deviation of background
208   Int_t count = 0;
209   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
210     if((injet[jpart] == -1) &&
211        (PtCutPass(jpart,nTracks))){
212       sigmaN += etbgTotalN/areaT - ptT[jpart];
213       // To be checked (Division by jet area to obtain standard deviation of rho ?)
214
215       count=count+1;
216     }
217   }
218   if (count>0)
219     sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
220
221 }
222
223 //----------------------------------------------------------------
224 void AliJetBkg::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN, Float_t&sigmaN,
225                                   const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
226                                   Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
227                                   Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet,
228                                   Int_t* injet, Float_t* &areaJet)
229 {
230   //
231   //background subtraction using statistical method
232   // Cases to take into account the EMCal geometry are included
233   //
234
235   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
236   Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
237   
238   //calculate energy inside
239   Float_t rc= header->GetRadius();
240   Float_t etIn[kMaxJets] = {0.0};
241   // Get number of tracks from EventCalTrk
242   Int_t nTracks = fEvent->GetNCalTrkTracks();
243   Float_t areaOut = 0.;
244
245   for(Int_t jpart = 0; jpart < nIn; jpart++)
246     { // loop for all particles in array
247       
248       for(Int_t ijet=0; ijet<nJ; ijet++)
249         {
250           Float_t deta = etaT[jpart] - etaJet[ijet];
251           Float_t dphi = phiT[jpart] - phiJet[ijet];
252           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
253           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
254           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
255           if(dr <= rc){ // particles inside this cone
256             if(jpart < nTracks) multJetT[ijet]++;
257             else multJetC[ijet]++;
258             multJet[ijet]++;
259             injet[jpart] = ijet;
260
261             if(PtCutPass(jpart,nTracks)){ // pt cut
262               etIn[ijet] += ptT[jpart];
263               if(SignalCutPass(jpart,nTracks))
264                 etsigJet[ijet]+= ptT[jpart];
265             }
266             break;
267           }
268         }// end jets loop
269     } //end particle loop
270   
271   // Calculate jet and background areas
272   Bool_t calcAreaOut = kFALSE;
273   CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut);
274
275   //subtract background using area method
276   for(Int_t ljet=0; ljet<nJ; ljet++){
277     Float_t areaRatio = areaJet[ljet]/areaOut;
278     etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
279   }
280   Int_t count=0;
281   etbgTotalN = etbgStat;
282
283   // estimate standard deviation of background
284   Float_t areaT = 0;
285   areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
286   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
287     if((injet[jpart] == -1) &&
288        (PtCutPass(jpart,nTracks))){
289       sigmaN += etbgTotalN/areaT - ptT[jpart];
290       count=count+1;
291     }
292   }
293   if(count>0)sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
294
295 }
296
297 //----------------------------------------------------------------
298 void AliJetBkg::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN,
299                                   const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, Float_t* etJet,
300                                   const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
301                                   Int_t* multJetT, Int_t* multJetC, Int_t* multJet, Int_t* injet, Float_t* &/*areaJet*/)
302 {
303   //
304   // Cone background subtraction method taking into acount dEt/deta distribution
305   // Cases to take into account the EMCal geometry are not included
306   //
307
308   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
309   //general
310   Float_t rc= header->GetRadius();
311   Float_t etamax = header->GetLegoEtaMax();
312   Float_t etamin = header->GetLegoEtaMin();
313   Int_t ndiv = 100;
314   // Get number of tracks from EventCalTrk
315   Int_t nTracks = fEvent->GetNCalTrkTracks();
316  
317   // jet energy and area arrays
318   Bool_t oldStatus = TH1::AddDirectoryStatus();
319   TH1::AddDirectory(kFALSE);
320
321   for(Int_t mjet=0; mjet<nJ; mjet++){
322     if(!fhEtJet[mjet]){
323       fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
324     }
325     if(!fhAreaJet[mjet]){
326       fhAreaJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
327     }
328     fhEtJet[mjet]->Reset();
329     fhAreaJet[mjet]->Reset();
330   }
331   // background energy and area
332   if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
333   fhEtBackg->Reset();
334   if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
335   fhAreaBackg->Reset();
336   TH1::AddDirectory(oldStatus);
337
338   //fill energies
339   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
340     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
341       Float_t deta = etaT[jpart] - etaJet[ijet];
342       Float_t dphi = phiT[jpart] - phiJet[ijet];
343       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
344       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
345       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
346       if(dr <= rc){ // particles inside this cone
347         if(jpart < nTracks) multJetT[ijet]++;
348         else multJetC[ijet]++;
349         multJet[ijet]++;
350         injet[jpart] = ijet;
351
352         if(PtCutPass(jpart,nTracks)){ // pt cut
353           fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
354           if(SignalCutPass(jpart,nTracks))
355             etsigJet[ijet]+= ptT[jpart];
356         }
357         break;
358       }
359     }// end jets loop
360
361     if((injet[jpart] == -1)  &&
362        (PtCutPass(jpart,nTracks) == 1))
363       fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
364   } //end particle loop
365
366   //calc areas
367   Float_t eta0 = etamin;
368   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
369   Float_t eta1 = eta0 + etaw;
370   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
371     Float_t etac = eta0 + etaw/2.0;
372     Float_t areabg = etaw*2.0*TMath::Pi();
373     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
374       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
375       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
376       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
377       Float_t areaj = 0.0;
378       if(deta0 > rc && deta1 < rc){
379         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
380         areaj = acc1;
381       }
382       if(deta0 < rc && deta1 > rc){
383         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
384         areaj = acc0;
385       }
386       if(deta0 < rc && deta1 < rc){
387         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
388         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
389         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
390         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
391         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
392       }
393       fhAreaJet[ijet]->Fill(etac,areaj);
394       areabg = areabg - areaj;
395     } // end jets loop
396     fhAreaBackg->Fill(etac,areabg);
397     eta0 = eta1;
398     eta1 = eta1 + etaw;
399   } // end loop for all eta bins
400
401   //subtract background
402   for(Int_t kjet=0; kjet<nJ; kjet++){
403     etJet[kjet] = 0.0; // first  clear etJet for this jet
404     for(Int_t bin = 0; bin< ndiv; bin++){
405       if(fhAreaJet[kjet]->GetBinContent(bin)){
406         Float_t areab = fhAreaBackg->GetBinContent(bin);
407         Float_t etb = fhEtBackg->GetBinContent(bin);
408         Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
409         etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
410       }
411     }
412   }
413
414   // calc background total
415   Double_t etOut = fhEtBackg->Integral();
416   Double_t areaOut = fhAreaBackg->Integral();
417   Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
418   etbgTotalN = etOut*areaT/areaOut;
419
420   Int_t count=0;
421   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
422     if((injet[jpart] == -1) &&
423        (PtCutPass(jpart,nTracks))){
424       sigmaN += etbgTotalN/areaT - ptT[jpart];
425       count=count+1;
426     }
427   }
428   sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
429   
430 }
431
432 //----------------------------------------------------------------
433 void AliJetBkg::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN,
434                                    const Float_t* ptT,const Float_t* etaT, const Float_t* phiT,
435                                    Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
436                                    Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet,
437                                    Int_t* injet,  Float_t* &/*areaJet*/)
438 {
439   // Ratio background subtraction method taking into acount dEt/deta distribution
440   // Cases to take into account the EMCal geometry are not included
441
442   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
443   //factor F calc before
444   Float_t bgRatioCut = header->GetBackgCutRatio();
445   
446   //general
447   Float_t rc= header->GetRadius();
448   Float_t etamax = header->GetLegoEtaMax();
449   Float_t etamin = header->GetLegoEtaMin();
450   Int_t ndiv = 100;
451   // Get number of tracks from EventCalTrk
452   Int_t nTracks = fEvent->GetNCalTrkTracks();
453   
454   // jet energy and area arrays
455   Bool_t oldStatus = TH1::AddDirectoryStatus();
456   TH1::AddDirectory(kFALSE);
457   for(Int_t mjet=0; mjet<nJ; mjet++){
458     if(!fhEtJet[mjet]){
459       fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
460     }
461     if(!fhAreaJet[mjet]){
462       fhAreaJet[mjet] = new TH1F(Form("hAreaJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
463     }
464     fhEtJet[mjet]->Reset();
465     fhAreaJet[mjet]->Reset();
466   }
467   // background energy and area
468   if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
469   fhEtBackg->Reset();
470   if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
471   fhAreaBackg->Reset();
472   TH1::AddDirectory(oldStatus);
473
474   //fill energies
475   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
476     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
477       Float_t deta = etaT[jpart] - etaJet[ijet];
478       Float_t dphi = phiT[jpart] - phiJet[ijet];
479       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
480       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
481       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
482       if(dr <= rc){ // particles inside this cone
483         if(jpart < nTracks) multJetT[ijet]++;
484         else multJetC[ijet]++;
485         multJet[ijet]++;
486         injet[jpart] = ijet;
487
488         if(PtCutPass(jpart,nTracks)){ // pt cut
489           fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
490           if(SignalCutPass(jpart,nTracks))
491             etsigJet[ijet]+= ptT[jpart];
492         }
493         break;
494       }
495     }// end jets loop
496     if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
497   } //end particle loop
498
499   //calc areas
500   Float_t eta0 = etamin;
501   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
502   Float_t eta1 = eta0 + etaw;
503   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
504     Float_t etac = eta0 + etaw/2.0;
505     Float_t areabg = etaw*2.0*TMath::Pi();
506     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
507       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
508       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
509       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
510       Float_t areaj = 0.0;
511       if(deta0 > rc && deta1 < rc){
512         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
513         areaj = acc1;
514       }
515       if(deta0 < rc && deta1 > rc){
516         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
517         areaj = acc0;
518       }
519       if(deta0 < rc && deta1 < rc){
520         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
521         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
522         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
523         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
524         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
525       }
526       fhAreaJet[ijet]->Fill(etac,areaj);
527       areabg = areabg - areaj;
528     } // end jets loop
529     fhAreaBackg->Fill(etac,areabg);
530     eta0 = eta1;
531     eta1 = eta1 + etaw;
532   } // end loop for all eta bins
533
534   //subtract background
535   for(Int_t kjet=0; kjet<nJ; kjet++){
536     etJet[kjet] = 0.0; // first  clear etJet for this jet
537     for(Int_t bin = 0; bin< ndiv; bin++){
538       if(fhAreaJet[kjet]->GetBinContent(bin)){
539         Float_t areab = fhAreaBackg->GetBinContent(bin);
540         Float_t etb = fhEtBackg->GetBinContent(bin);
541         Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
542         etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
543       }
544     }
545   }
546
547   // calc background total
548   Double_t etOut = fhEtBackg->Integral();
549   Double_t areaOut = fhAreaBackg->Integral();
550   Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
551   etbgTotalN = etOut*areaT/areaOut;
552  
553   Int_t count=0;
554     
555   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
556     if((injet[jpart] == -1) &&
557        (PtCutPass(jpart,nTracks))){
558       sigmaN += etbgTotalN/areaT - ptT[jpart];
559       count=count+1;
560     }
561   }
562   sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
563
564 }
565