JETAN/DEV/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, const 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, const 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   sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
 219
 220 }
 221
 222 //----------------------------------------------------------------
 223 void AliJetBkg::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN, Float_t&sigmaN,
 224                                   const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
 225                                   Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
 226                                   Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet,
 227                                   Int_t* injet, Float_t* &areaJet)
 228 {
 229   //
 230   //background subtraction using statistical method
 231   // Cases to take into account the EMCal geometry are included
 232   //
 233
 234   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
 235   Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
 236
 237   //calculate energy inside
 238   Float_t rc= header->GetRadius();
 239   Float_t etIn[kMaxJets] = {0.0};
 240   // Get number of tracks from EventCalTrk
 241   Int_t nTracks = fEvent->GetNCalTrkTracks();
 242   Float_t areaOut = 0.;
 243
 244   for(Int_t jpart = 0; jpart < nIn; jpart++)
 245     { // loop for all particles in array
 246
 247       for(Int_t ijet=0; ijet<nJ; ijet++)
 248         {
 249           Float_t deta = etaT[jpart] - etaJet[ijet];
 250           Float_t dphi = phiT[jpart] - phiJet[ijet];
 251           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 252           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 253           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 254           if(dr <= rc){ // particles inside this cone
 255             if(jpart < nTracks) multJetT[ijet]++;
 256             else multJetC[ijet]++;
 257             multJet[ijet]++;
 258             injet[jpart] = ijet;
 259
 260             if(PtCutPass(jpart,nTracks)){ // pt cut
 261               etIn[ijet] += ptT[jpart];
 262               if(SignalCutPass(jpart,nTracks))
 263                 etsigJet[ijet]+= ptT[jpart];
 264             }
 265             break;
 266           }
 267         }// end jets loop
 268     } //end particle loop
 269
 270   // Calculate jet and background areas
 271   Bool_t calcAreaOut = kFALSE;
 272   CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut);
 273
 274   //subtract background using area method
 275   for(Int_t ljet=0; ljet<nJ; ljet++){
 276     Float_t areaRatio = areaJet[ljet]/areaOut;
 277     etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
 278   }
 279   Int_t count=0;
 280   etbgTotalN = etbgStat;
 281
 282   // estimate standard deviation of background
 283   Float_t areaT = 0;
 284   areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
 285   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 286     if((injet[jpart] == -1) &&
 287        (PtCutPass(jpart,nTracks))){
 288       sigmaN += etbgTotalN/areaT - ptT[jpart];
 289       count=count+1;
 290     }
 291   }
 292   sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
 293
 294 }
 295
 296 //----------------------------------------------------------------
 297 void AliJetBkg::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN,
 298                                   const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, Float_t* etJet,
 299                                   const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
 300                                   Int_t* multJetT, Int_t* multJetC, Int_t* multJet, Int_t* injet, Float_t* &/*areaJet*/)
 301 {
 302   //
 303   // Cone background subtraction method taking into acount dEt/deta distribution
 304   // Cases to take into account the EMCal geometry are not included
 305   //
 306
 307   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
 308   //general
 309   Float_t rc= header->GetRadius();
 310   Float_t etamax = header->GetLegoEtaMax();
 311   Float_t etamin = header->GetLegoEtaMin();
 312   Int_t ndiv = 100;
 313   // Get number of tracks from EventCalTrk
 314   Int_t nTracks = fEvent->GetNCalTrkTracks();
 315
 316   // jet energy and area arrays
 317   Bool_t oldStatus = TH1::AddDirectoryStatus();
 318   TH1::AddDirectory(kFALSE);
 319
 320   for(Int_t mjet=0; mjet<nJ; mjet++){
 321     if(!fhEtJet[mjet]){
 322       fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
 323     }
 324     if(!fhAreaJet[mjet]){
 325       fhAreaJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
 326     }
 327     fhEtJet[mjet]->Reset();
 328     fhAreaJet[mjet]->Reset();
 329   }
 330   // background energy and area
 331   if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
 332   fhEtBackg->Reset();
 333   if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
 334   fhAreaBackg->Reset();
 335   TH1::AddDirectory(oldStatus);
 336
 337   //fill energies
 338   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 339     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 340       Float_t deta = etaT[jpart] - etaJet[ijet];
 341       Float_t dphi = phiT[jpart] - phiJet[ijet];
 342       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 343       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 344       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 345       if(dr <= rc){ // particles inside this cone
 346         if(jpart < nTracks) multJetT[ijet]++;
 347         else multJetC[ijet]++;
 348         multJet[ijet]++;
 349         injet[jpart] = ijet;
 350
 351         if(PtCutPass(jpart,nTracks)){ // pt cut
 352           fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
 353           if(SignalCutPass(jpart,nTracks))
 354             etsigJet[ijet]+= ptT[jpart];
 355         }
 356         break;
 357       }
 358     }// end jets loop
 359
 360     if((injet[jpart] == -1)  &&
 361        (PtCutPass(jpart,nTracks) == 1))
 362       fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
 363   } //end particle loop
 364
 365   //calc areas
 366   Float_t eta0 = etamin;
 367   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
 368   Float_t eta1 = eta0 + etaw;
 369   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
 370     Float_t etac = eta0 + etaw/2.0;
 371     Float_t areabg = etaw*2.0*TMath::Pi();
 372     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 373       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
 374       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
 375       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
 376       Float_t areaj = 0.0;
 377       if(deta0 > rc && deta1 < rc){
 378         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 379         areaj = acc1;
 380       }
 381       if(deta0 < rc && deta1 > rc){
 382         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 383         areaj = acc0;
 384       }
 385       if(deta0 < rc && deta1 < rc){
 386         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 387         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 388         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
 389         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
 390         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
 391       }
 392       fhAreaJet[ijet]->Fill(etac,areaj);
 393       areabg = areabg - areaj;
 394     } // end jets loop
 395     fhAreaBackg->Fill(etac,areabg);
 396     eta0 = eta1;
 397     eta1 = eta1 + etaw;
 398   } // end loop for all eta bins
 399
 400   //subtract background
 401   for(Int_t kjet=0; kjet<nJ; kjet++){
 402     etJet[kjet] = 0.0; // first  clear etJet for this jet
 403     for(Int_t bin = 0; bin< ndiv; bin++){
 404       if(fhAreaJet[kjet]->GetBinContent(bin)){
 405         Float_t areab = fhAreaBackg->GetBinContent(bin);
 406         Float_t etb = fhEtBackg->GetBinContent(bin);
 407         Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
 408         etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
 409       }
 410     }
 411   }
 412
 413   // calc background total
 414   Double_t etOut = fhEtBackg->Integral();
 415   Double_t areaOut = fhAreaBackg->Integral();
 416   Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
 417   etbgTotalN = etOut*areaT/areaOut;
 418
 419   Int_t count=0;
 420   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 421     if((injet[jpart] == -1) &&
 422        (PtCutPass(jpart,nTracks))){
 423       sigmaN += etbgTotalN/areaT - ptT[jpart];
 424       count=count+1;
 425     }
 426   }
 427   sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
 428
 429 }
 430
 431 //----------------------------------------------------------------
 432 void AliJetBkg::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN,
 433                                    const Float_t* ptT,const Float_t* etaT, const Float_t* phiT,
 434                                    Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
 435                                    Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet,
 436                                    Int_t* injet,  Float_t* &/*areaJet*/)
 437 {
 438   // Ratio background subtraction method taking into acount dEt/deta distribution
 439   // Cases to take into account the EMCal geometry are not included
 440
 441   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
 442   //factor F calc before
 443   Float_t bgRatioCut = header->GetBackgCutRatio();
 444
 445   //general
 446   Float_t rc= header->GetRadius();
 447   Float_t etamax = header->GetLegoEtaMax();
 448   Float_t etamin = header->GetLegoEtaMin();
 449   Int_t ndiv = 100;
 450   // Get number of tracks from EventCalTrk
 451   Int_t nTracks = fEvent->GetNCalTrkTracks();
 452
 453   // jet energy and area arrays
 454   Bool_t oldStatus = TH1::AddDirectoryStatus();
 455   TH1::AddDirectory(kFALSE);
 456   for(Int_t mjet=0; mjet<nJ; mjet++){
 457     if(!fhEtJet[mjet]){
 458       fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
 459     }
 460     if(!fhAreaJet[mjet]){
 461       fhAreaJet[mjet] = new TH1F(Form("hAreaJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
 462     }
 463     fhEtJet[mjet]->Reset();
 464     fhAreaJet[mjet]->Reset();
 465   }
 466   // background energy and area
 467   if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
 468   fhEtBackg->Reset();
 469   if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
 470   fhAreaBackg->Reset();
 471   TH1::AddDirectory(oldStatus);
 472
 473   //fill energies
 474   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 475     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 476       Float_t deta = etaT[jpart] - etaJet[ijet];
 477       Float_t dphi = phiT[jpart] - phiJet[ijet];
 478       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 479       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 480       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 481       if(dr <= rc){ // particles inside this cone
 482         if(jpart < nTracks) multJetT[ijet]++;
 483         else multJetC[ijet]++;
 484         multJet[ijet]++;
 485         injet[jpart] = ijet;
 486
 487         if(PtCutPass(jpart,nTracks)){ // pt cut
 488           fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
 489           if(SignalCutPass(jpart,nTracks))
 490             etsigJet[ijet]+= ptT[jpart];
 491         }
 492         break;
 493       }
 494     }// end jets loop
 495     if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
 496   } //end particle loop
 497
 498   //calc areas
 499   Float_t eta0 = etamin;
 500   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
 501   Float_t eta1 = eta0 + etaw;
 502   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
 503     Float_t etac = eta0 + etaw/2.0;
 504     Float_t areabg = etaw*2.0*TMath::Pi();
 505     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 506       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
 507       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
 508       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
 509       Float_t areaj = 0.0;
 510       if(deta0 > rc && deta1 < rc){
 511         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 512         areaj = acc1;
 513       }
 514       if(deta0 < rc && deta1 > rc){
 515         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 516         areaj = acc0;
 517       }
 518       if(deta0 < rc && deta1 < rc){
 519         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 520         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 521         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
 522         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
 523         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
 524       }
 525       fhAreaJet[ijet]->Fill(etac,areaj);
 526       areabg = areabg - areaj;
 527     } // end jets loop
 528     fhAreaBackg->Fill(etac,areabg);
 529     eta0 = eta1;
 530     eta1 = eta1 + etaw;
 531   } // end loop for all eta bins
 532
 533   //subtract background
 534   for(Int_t kjet=0; kjet<nJ; kjet++){
 535     etJet[kjet] = 0.0; // first  clear etJet for this jet
 536     for(Int_t bin = 0; bin< ndiv; bin++){
 537       if(fhAreaJet[kjet]->GetBinContent(bin)){
 538         Float_t areab = fhAreaBackg->GetBinContent(bin);
 539         Float_t etb = fhEtBackg->GetBinContent(bin);
 540         Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
 541         etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
 542       }
 543     }
 544   }
 545
 546   // calc background total
 547   Double_t etOut = fhEtBackg->Integral();
 548   Double_t areaOut = fhAreaBackg->Integral();
 549   Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin());
 550   etbgTotalN = etOut*areaT/areaOut;
 551
 552   Int_t count=0;
 553
 554   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 555     if((injet[jpart] == -1) &&
 556        (PtCutPass(jpart,nTracks))){
 557       sigmaN += etbgTotalN/areaT - ptT[jpart];
 558       count=count+1;
 559     }
 560   }
 561   sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count);
 562
 563 }
 564