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