+ AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
+ //factor F calc before
+ Float_t bgRatioCut = header->GetBackgCutRatio();
+
+ //general
+ Float_t rc= header->GetRadius();
+ Float_t etamax = header->GetLegoEtaMax();
+ Float_t etamin = header->GetLegoEtaMin();
+ Int_t ndiv = 100;
+ // Get number of tracks from EventCalTrk
+ Int_t nTracks = fEvent->GetNCalTrkTracks();
+
+ // jet energy and area arrays
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+ for(Int_t mjet=0; mjet<nJ; mjet++){
+ if(!fhEtJet[mjet]){
+ fhEtJet[mjet] = new TH1F(Form("hEtJet%d", mjet),"et dist in eta ",ndiv,etamin,etamax);
+ }
+ if(!fhAreaJet[mjet]){
+ fhAreaJet[mjet] = new TH1F(Form("hAreaJet%d", mjet),"area dist in eta ",ndiv,etamin,etamax);
+ }
+ fhEtJet[mjet]->Reset();
+ fhAreaJet[mjet]->Reset();
+ }
+ // background energy and area
+ if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
+ fhEtBackg->Reset();
+ if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
+ fhAreaBackg->Reset();
+ TH1::AddDirectory(oldStatus);
+
+ //fill energies
+ for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+ for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
+ Float_t deta = etaT[jpart] - etaJet[ijet];
+ Float_t dphi = phiT[jpart] - phiJet[ijet];
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if(dr <= rc){ // particles inside this cone
+ if(jpart < nTracks) multJetT[ijet]++;
+ else multJetC[ijet]++;
+ multJet[ijet]++;
+ injet[jpart] = ijet;
+
+ if(PtCutPass(jpart,nTracks)){ // pt cut
+ fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
+ if(SignalCutPass(jpart,nTracks))
+ etsigJet[ijet]+= ptT[jpart];
+ }
+ break;
+ }
+ }// end jets loop
+ if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+ } //end particle loop
+
+ //calc areas
+ Float_t eta0 = etamin;
+ Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
+ Float_t eta1 = eta0 + etaw;
+ for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
+ Float_t etac = eta0 + etaw/2.0;
+ Float_t areabg = etaw*2.0*TMath::Pi();
+ for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
+ Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
+ Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
+ Float_t acc0 = 0.0; Float_t acc1 = 0.0;
+ Float_t areaj = 0.0;
+ if(deta0 > rc && deta1 < rc){
+ acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+ areaj = acc1;
+ }
+ if(deta0 < rc && deta1 > rc){
+ acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+ areaj = acc0;
+ }
+ if(deta0 < rc && deta1 < rc){
+ acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+ acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+ if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
+ if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
+ if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
+ }
+ fhAreaJet[ijet]->Fill(etac,areaj);
+ areabg = areabg - areaj;
+ } // end jets loop
+ fhAreaBackg->Fill(etac,areabg);
+ eta0 = eta1;
+ eta1 = eta1 + etaw;
+ } // end loop for all eta bins
+
+ //subtract background
+ for(Int_t kjet=0; kjet<nJ; kjet++){
+ etJet[kjet] = 0.0; // first clear etJet for this jet
+ for(Int_t bin = 0; bin< ndiv; bin++){
+ if(fhAreaJet[kjet]->GetBinContent(bin)){
+ Float_t areab = fhAreaBackg->GetBinContent(bin);
+ Float_t etb = fhEtBackg->GetBinContent(bin);
+ Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab;
+ etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
+ }
+ }
+ }