1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //---------------------------------------------------------------------
19 // UA1 Cone Algorithm Jet finder for charged + neutral jet studies
20 // manages the search for jets using charged particle momentum and
21 // neutral cell energy information
22 // Based on UA1 V1 (from R. Diaz)
23 // Author: magali.estienne@subatech.in2p3.fr
24 //---------------------------------------------------------------------
26 #include <TClonesArray.h>
29 #include <TLorentzVector.h>
31 #include <TRefArray.h>
34 #include "AliUA1JetFinderV2.h"
35 #include "AliUA1JetHeaderV1.h"
36 #include "AliJetUnitArray.h"
43 ClassImp(AliUA1JetFinderV2)
46 ////////////////////////////////////////////////////////////////////////
47 AliUA1JetFinderV2::AliUA1JetFinderV2() :
58 ////////////////////////////////////////////////////////////////////////
59 AliUA1JetFinderV2::~AliUA1JetFinderV2()
66 ////////////////////////////////////////////////////////////////////////
67 void AliUA1JetFinderV2::FindJetsC()
70 // Used to find jets using charged particle momentum information
72 // 1) Fill cell map array
73 // 2) calculate total energy and fluctuation level
75 // 3.1) look centroides in cell map
76 // 3.2) calculate total energy in cones
77 // 3.3) flag as a possible jet
78 // 3.4) reorder cones by energy
79 // 4) subtract backg in accepted jets
80 // 5) fill AliJet list
82 // Transform input to pt,eta,phi plus lego
84 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
85 TClonesArray* lvArray = fReader->GetMomentumArray();
86 Int_t nIn = lvArray->GetEntries();
90 // local arrays for input
91 Float_t* ptT = new Float_t[nIn];
92 Float_t* etaT = new Float_t[nIn];
93 Float_t* phiT = new Float_t[nIn];
94 Float_t* cFlagT = new Float_t[nIn]; // Temporarily added
95 Float_t* sFlagT = new Float_t[nIn]; // Temporarily added
96 Int_t* injet = new Int_t[nIn];
98 //total energy in array
99 Float_t etbgTotal = 0.0;
100 TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
102 // load input vectors and calculate total energy in array
103 for (Int_t i = 0; i < nIn; i++){
104 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
107 phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
108 cFlagT[i] = fReader->GetCutFlag(i);
109 sFlagT[i] = fReader->GetSignalFlag(i);
111 if (fReader->GetCutFlag(i) != 1) continue;
112 fLego->Fill(etaT[i], phiT[i], ptT[i]);
113 hPtTotal->Fill(ptT[i]);
117 fJets->SetNinput(nIn);
119 // calculate total energy and fluctuation in map
120 Double_t meanpt = hPtTotal->GetMean();
121 Double_t ptRMS = hPtTotal->GetRMS();
122 Double_t npart = hPtTotal->GetEntries();
123 Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
125 // arrays to hold jets
126 Float_t* etaJet = new Float_t[30]; // eta jet
127 Float_t* phiJet = new Float_t[30]; // phi jet
128 Float_t* etJet = new Float_t[30]; // et jet
129 Float_t* etsigJet = new Float_t[30]; // signal et in jet
130 Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable)
131 Int_t* ncellsJet = new Int_t[30];
132 Int_t* multJet = new Int_t[30];
133 //--- Added for jet reordering at the end of the jet finding procedure
134 Float_t* etaJetOk = new Float_t[30];
135 Float_t* phiJetOk = new Float_t[30];
136 Float_t* etJetOk = new Float_t[30];
137 Float_t* etsigJetOk = new Float_t[30]; // signal et in jet
138 Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable)
139 Int_t* ncellsJetOk = new Int_t[30];
140 Int_t* multJetOk = new Int_t[30];
141 //--------------------------
142 Int_t nJets; // to hold number of jets found by algorithm
143 Int_t nj; // number of jets accepted
144 Float_t prec = header->GetPrecBg();
146 while(bgprec > prec){
147 //reset jet arrays in memory
148 memset(etaJet,0,sizeof(Float_t)*30);
149 memset(phiJet,0,sizeof(Float_t)*30);
150 memset(etJet,0,sizeof(Float_t)*30);
151 memset(etallJet,0,sizeof(Float_t)*30);
152 memset(etsigJet,0,sizeof(Float_t)*30);
153 memset(ncellsJet,0,sizeof(Int_t)*30);
154 memset(multJet,0,sizeof(Int_t)*30);
155 //--- Added for jet reordering at the end of the jet finding procedure
156 memset(etaJetOk,0,sizeof(Float_t)*30);
157 memset(phiJetOk,0,sizeof(Float_t)*30);
158 memset(etJetOk,0,sizeof(Float_t)*30);
159 memset(etallJetOk,0,sizeof(Float_t)*30);
160 memset(etsigJetOk,0,sizeof(Float_t)*30);
161 memset(ncellsJetOk,0,sizeof(Int_t)*30);
162 memset(multJetOk,0,sizeof(Int_t)*30);
163 //--------------------------
167 // reset particles-jet array in memory
168 memset(injet,-1,sizeof(Int_t)*nIn);
169 //run cone algorithm finder
170 RunAlgoritmC(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
172 //run background subtraction
173 if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
174 nj = header->GetNAcceptJets();
177 //subtract background
178 Float_t etbgTotalN = 0.0; //new background
179 if(header->GetBackgMode() == 1) // standard
180 SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
181 if(header->GetBackgMode() == 2) //cone
182 SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
183 if(header->GetBackgMode() == 3) //ratio
184 SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
185 if(header->GetBackgMode() == 4) //statistic
186 SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
188 if(etbgTotalN != 0.0)
189 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
192 etbgTotal = etbgTotalN; // update with new background estimation
196 Int_t* idxjets = new Int_t[nj];
198 printf("Found %d jets \n", nj);
200 // Reorder jets by et in cone
201 Int_t * idx = new Int_t[nJets];
202 TMath::Sort(nJets, etJet, idx);
203 for(Int_t p = 0; p < nJets; p++){
204 etaJetOk[p] = etaJet[idx[p]];
205 phiJetOk[p] = phiJet[idx[p]];
206 etJetOk[p] = etJet[idx[p]];
207 etallJetOk[p] = etJet[idx[p]];
208 etsigJetOk[p] = etsigJet[idx[p]];
209 ncellsJetOk[p] = ncellsJet[idx[p]];
210 multJetOk[p] = multJet[idx[p]];
213 for(Int_t kj=0; kj<nj; kj++)
215 if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
216 (etaJetOk[kj] < (header->GetJetEtaMin())) ||
217 (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
218 Float_t px, py,pz,en; // convert to 4-vector
219 px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
220 py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
221 pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
222 en = TMath::Sqrt(px * px + py * py + pz * pz);
223 fJets->AddJet(px, py, pz, en);
225 AliAODJet jet(px, py, pz, en);
230 idxjets[nselectj] = kj;
234 //add signal percentage and total signal in AliJets for analysis tool
235 Float_t* percentage = new Float_t[nselectj];
236 Int_t* ncells = new Int_t[nselectj];
237 Int_t* mult = new Int_t[nselectj];
238 for(Int_t i = 0; i< nselectj; i++)
240 percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
241 ncells[i] = ncellsJetOk[idxjets[i]];
242 mult[i] = multJetOk[idxjets[i]];
245 //add particle-injet relationship ///
246 for(Int_t bj = 0; bj < nIn; bj++)
248 if(injet[bj] == -1) continue; //background particle
250 for(Int_t ci = 0; ci< nselectj; ci++){
251 if(injet[bj] == idxjets[ci]){
257 if(bflag == 0) injet[bj] = -1; // set as background particle
260 fJets->SetNCells(ncells);
261 fJets->SetPtFromSignal(percentage);
262 fJets->SetMultiplicities(mult);
263 fJets->SetInJet(injet);
264 fJets->SetEtaIn(etaT);
265 fJets->SetPhiIn(phiT);
267 fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi()));
288 //--- Added for jet reordering
296 //--------------------------
300 ////////////////////////////////////////////////////////////////////////
301 void AliUA1JetFinderV2::FindJets()
304 // Used to find jets using charged particle momentum information
305 // & neutral energy from calo cells
307 // 1) Fill cell map array
308 // 2) calculate total energy and fluctuation level
310 // 3.1) look centroides in cell map
311 // 3.2) calculate total energy in cones
312 // 3.3) flag as a possible jet
313 // 3.4) reorder cones by energy
314 // 4) subtract backg in accepted jets
315 // 5) fill AliJet list
317 // transform input to pt,eta,phi plus lego
319 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
320 TClonesArray* fUnit = fReader->GetUnitArray();
321 Int_t nCand = fReader->GetNumCandidate();
322 Int_t nCandCut = fReader->GetNumCandidateCut();
323 Int_t nIn = fUnit->GetEntries();
324 Float_t ptMin = fReader->GetReaderHeader()->GetPtCut();
326 fDebug = fReader->GetReaderHeader()->GetDebug();
327 if (nIn == 0) return;
329 Int_t nCandidateCut = 0;
330 Int_t nCandidate = 0;
333 nCandidateCut = nCandCut;
335 // local arrays for input No Cuts
336 // Both pt < ptMin and pt > ptMin
337 Float_t* ptT = new Float_t[nCandidate];
338 Float_t* en2T = new Float_t[nCandidate];
339 Float_t* pt2T = new Float_t[nCandidate];
340 Int_t* detT = new Int_t[nCandidate];
341 Float_t* etaT = new Float_t[nCandidate];
342 Float_t* phiT = new Float_t[nCandidate];
343 Float_t* cFlagT = new Float_t[nCandidate];
344 Float_t* cFlag2T = new Float_t[nCandidate];
345 Float_t* sFlagT = new Float_t[nCandidate];
346 Float_t* cClusterT = new Float_t[nCandidate];
347 Int_t* vectT = new Int_t[nCandidate];
349 Int_t* injet = new Int_t[nCandidate];
350 Int_t* sflag = new Int_t[nCandidate];
351 TRefArray* trackRef = new TRefArray();
353 //total energy in array
354 Float_t etbgTotal = 0.0;
355 TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
358 Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
359 Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
360 Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
361 Int_t *flagCell = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray
362 Float_t *etCell2 = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
363 Float_t *etaCell2 = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
364 Float_t *phiCell2 = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
365 Int_t *flagCell2 = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray
367 // Information extracted from fUnitArray
368 // Load input vectors and calculate total energy in array
369 for(Int_t i=0; i<nIn; i++)
371 // Recover particle information from UnitArray
373 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
374 TRefArray* ref = uArray->GetUnitTrackRef();
375 Int_t nRef = ref->GetEntries();
377 if(uArray->GetUnitEnergy()>0.){
379 for(Int_t jpart=0; jpart<nRef;jpart++)
380 trackRef->Add((AliVTrack*)ref->At(jpart));
381 ptT[loop1] = uArray->GetUnitEnergy();
382 detT[loop1] = uArray->GetUnitDetectorFlag();
383 etaT[loop1] = uArray->GetUnitEta();
384 phiT[loop1] = uArray->GetUnitPhi();
385 cFlagT[loop1]= uArray->GetUnitCutFlag(); // pt cut tpc
386 cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal
387 sFlagT[loop1]= uArray->GetUnitSignalFlag();
389 if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) {
393 en2T[loop1] = ptT[loop1] - header->GetMinCellEt();
394 if(en2T[loop1] < 0) en2T[loop1]=0;
395 hPtTotal->Fill(en2T[loop1]);
396 etbgTotal += en2T[loop1];
398 if(detT[loop1]==0){ // TPC+ITS
400 for(Int_t j=0; j<nRef;j++){
401 Float_t x=0.; Float_t y=0.; Float_t z=0.;
402 x = ((AliVTrack*)ref->At(j))->Px();
403 y = ((AliVTrack*)ref->At(j))->Py();
404 z = ((AliVTrack*)ref->At(j))->Pz();
405 pt = TMath::Sqrt(x*x+y*y);
414 if(detT[loop1]==2) { // EMCal
418 for(Int_t j=0; j<nRef;j++){
419 Float_t x=0.; Float_t y=0.; Float_t z=0.;
420 x = ((AliVTrack*)ref->At(j))->Px();
421 y = ((AliVTrack*)ref->At(j))->Py();
422 z = ((AliVTrack*)ref->At(j))->Pz();
423 pt = TMath::Sqrt(x*x+y*y);
432 enC = ptT[loop1] - ptCTot - header->GetMinCellEt();
442 if(uArray->GetUnitCutFlag()==1) {
443 if(uArray->GetUnitDetectorFlag()==1){ // EMCal case
444 etCell[i] = uArray->GetUnitEnergy() - header->GetMinCellEt();
445 if ((uArray->GetUnitEnergy() - header->GetMinCellEt()) < 0.0) etCell[i]=0.;
446 etaCell[i] = uArray->GetUnitEta();
447 phiCell[i] = uArray->GetUnitPhi();
448 flagCell[i] = 0; // default
449 etCell2[i] = etCell[i];
450 etaCell2[i] = uArray->GetUnitEta();
451 phiCell2[i] = uArray->GetUnitPhi();
452 flagCell2[i] = 0; // default
454 if(uArray->GetUnitDetectorFlag()==0){ // TPC case
455 Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
456 for(Int_t j=0; j<nRef;j++)
458 Float_t x=0.; Float_t y=0.; Float_t z=0.;
459 x = ((AliVTrack*)ref->At(j))->Px();
460 y = ((AliVTrack*)ref->At(j))->Py();
461 z = ((AliVTrack*)ref->At(j))->Pz();
462 pt = TMath::Sqrt(x*x+y*y);
470 if(et1 < 0.) etCell[i] = etCell2[i] = 0.;
471 etaCell[i] = uArray->GetUnitEta();
472 phiCell[i] = uArray->GetUnitPhi();
473 flagCell[i] = 0; // default
474 etaCell2[i] = uArray->GetUnitEta();
475 phiCell2[i] = uArray->GetUnitPhi();
476 flagCell2[i] = 0; // default
478 if(uArray->GetUnitDetectorFlag()==2){ // TPC + EMCal case
480 Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
482 for(Int_t j=0; j<nRef;j++)
484 Float_t x=0.; Float_t y=0.; Float_t z=0.;
485 x = ((AliVTrack*)ref->At(j))->Px();
486 y = ((AliVTrack*)ref->At(j))->Py();
487 z = ((AliVTrack*)ref->At(j))->Pz();
488 pt = TMath::Sqrt(x*x+y*y);
495 enC = uArray->GetUnitEnergy() - ptCTot;
496 etCell[i] = et1 + enC - header->GetMinCellEt();
497 etCell2[i] = et2 + enC - header->GetMinCellEt();
498 if((enC + et1 - header->GetMinCellEt()) < 0.) etCell[i] = etCell2[i] = 0.;
499 etaCell[i] = uArray->GetUnitEta();
500 phiCell[i] = uArray->GetUnitPhi();
501 flagCell[i] = 0; // default
502 etaCell2[i] = uArray->GetUnitEta();
503 phiCell2[i] = uArray->GetUnitPhi();
504 flagCell2[i] = 0; // default
509 etaCell[i] = uArray->GetUnitEta();
510 phiCell[i] = uArray->GetUnitPhi();
513 etaCell2[i] = uArray->GetUnitEta();
514 phiCell2[i] = uArray->GetUnitPhi();
517 } // end loop on nCandidate
519 fJets->SetNinput(nCandidate);
521 // calculate total energy and fluctuation in map
522 Double_t meanpt = hPtTotal->GetMean();
523 Double_t ptRMS = hPtTotal->GetRMS();
524 Double_t npart = hPtTotal->GetEntries();
525 Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
527 // arrays to hold jets
528 Float_t* etaJet = new Float_t[30];
529 Float_t* phiJet = new Float_t[30];
530 Float_t* etJet = new Float_t[30];
531 Float_t* etsigJet = new Float_t[30]; //signal et in jet
532 Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable)
533 Int_t* ncellsJet = new Int_t[30];
534 Int_t* multJet = new Int_t[30];
535 //--- Added by me for jet reordering at the end of the jet finding procedure
536 Float_t* etaJetOk = new Float_t[30];
537 Float_t* phiJetOk = new Float_t[30];
538 Float_t* etJetOk = new Float_t[30];
539 Float_t* etsigJetOk = new Float_t[30]; //signal et in jet
540 Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable)
541 Int_t* ncellsJetOk = new Int_t[30];
542 Int_t* multJetOk = new Int_t[30];
543 //--------------------------
544 Int_t nJets; // to hold number of jets found by algorithm
545 Int_t nj; // number of jets accepted
546 Float_t prec = header->GetPrecBg();
549 while(bgprec > prec){
551 //reset jet arrays in memory
552 memset(etaJet,0,sizeof(Float_t)*30);
553 memset(phiJet,0,sizeof(Float_t)*30);
554 memset(etJet,0,sizeof(Float_t)*30);
555 memset(etallJet,0,sizeof(Float_t)*30);
556 memset(etsigJet,0,sizeof(Float_t)*30);
557 memset(ncellsJet,0,sizeof(Int_t)*30);
558 memset(multJet,0,sizeof(Int_t)*30);
559 //--- Added by me for jet reordering at the end of the jet finding procedure
560 memset(etaJetOk,0,sizeof(Float_t)*30);
561 memset(phiJetOk,0,sizeof(Float_t)*30);
562 memset(etJetOk,0,sizeof(Float_t)*30);
563 memset(etallJetOk,0,sizeof(Float_t)*30);
564 memset(etsigJetOk,0,sizeof(Float_t)*30);
565 memset(ncellsJetOk,0,sizeof(Int_t)*30);
566 memset(multJetOk,0,sizeof(Int_t)*30);
571 // reset particles-jet array in memory
572 memset(injet,-1,sizeof(Int_t)*nCandidate);
573 //run cone algorithm finder
574 RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etCell2,etaCell2,phiCell2,
575 flagCell2,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,
578 //run background subtraction
579 if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
580 nj = header->GetNAcceptJets();
584 //subtract background
585 Float_t etbgTotalN = 0.0; //new background
586 if(header->GetBackgMode() == 1) // standard
587 SubtractBackg(nCandidate,nj,etbgTotalN,en2T,vectT,etaT,phiT,cFlagT,cFlag2T,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
588 // To be modified ------------------------
589 if(header->GetBackgMode() == 2) //cone
590 SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
591 if(header->GetBackgMode() == 3) //ratio
592 SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
593 if(header->GetBackgMode() == 4) //statistic
594 SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
595 //----------------------------------------
597 if(etbgTotalN != 0.0)
598 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
601 etbgTotal = etbgTotalN; // update with new background estimation
605 Int_t* idxjets = new Int_t[nj];
607 printf("Found %d jets \n", nj);
609 // Reorder jets by et in cone
610 // Sort jets by energy
611 Int_t * idx = new Int_t[nJets];
612 TMath::Sort(nJets, etJet, idx);
613 for(Int_t p = 0; p < nJets; p++)
615 etaJetOk[p] = etaJet[idx[p]];
616 phiJetOk[p] = phiJet[idx[p]];
617 etJetOk[p] = etJet[idx[p]];
618 etallJetOk[p] = etJet[idx[p]];
619 etsigJetOk[p] = etsigJet[idx[p]];
620 ncellsJetOk[p] = ncellsJet[idx[p]];
621 multJetOk[p] = multJet[idx[p]];
624 for(Int_t kj=0; kj<nj; kj++)
626 if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
627 (etaJetOk[kj] < (header->GetJetEtaMin())) ||
628 (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
629 Float_t px, py,pz,en; // convert to 4-vector
630 px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
631 py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
632 pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
633 en = TMath::Sqrt(px * px + py * py + pz * pz);
634 fJets->AddJet(px, py, pz, en);
635 AliAODJet jet(px, py, pz, en);
640 idxjets[nselectj] = kj;
644 //add signal percentage and total signal in AliJets for analysis tool
645 Float_t* percentage = new Float_t[nselectj];
646 Int_t* ncells = new Int_t[nselectj];
647 Int_t* mult = new Int_t[nselectj];
648 for(Int_t i = 0; i< nselectj; i++)
650 percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
651 ncells[i] = ncellsJetOk[idxjets[i]];
652 mult[i] = multJetOk[idxjets[i]];
655 //add particle-injet relationship ///
656 for(Int_t bj = 0; bj < nCandidate; bj++)
658 if(injet[bj] == -1) continue; //background particle
660 for(Int_t ci = 0; ci< nselectj; ci++){
661 if(injet[bj] == idxjets[ci]){
667 if(bflag == 0) injet[bj] = -1; // set as background particle
670 fJets->SetNCells(ncells);
671 fJets->SetPtFromSignal(percentage);
672 fJets->SetMultiplicities(mult);
673 fJets->SetInJet(injet);
674 fJets->SetEtaIn(etaT);
675 fJets->SetPhiIn(phiT);
677 fJets->SetTrackReferences(trackRef);
678 fJets->SetDetectorFlagIn(detT);
679 fJets->SetEtAvg(etbgTotal/(2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax()-header->GetLegoPhiMin())));
713 //--- Added for jet reordering
721 //--------------------------
729 ////////////////////////////////////////////////////////////////////////
730 void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell,
731 Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2,
732 const Float_t* phiCell2, const Int_t* flagCell2, Float_t etbgTotal,
733 Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet,
734 Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet)
737 // Main method for jet finding
738 // UA1 base cone finder
742 fDebug = fReader->GetReaderHeader()->GetDebug();
745 // Check enough space! *to be done*
746 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
747 for(Int_t i=0; i<nCell; i++){
748 etCell[i] = etCell2[i];
749 etaCell[i] = etaCell2[i];
750 phiCell[i] = phiCell2[i];
751 flagCell[i] = flagCell2[i];
754 // Parameters from header
755 Float_t minmove = header->GetMinMove();
756 Float_t maxmove = header->GetMaxMove();
757 Float_t rc = header->GetRadius();
758 Float_t etseed = header->GetEtSeed();
760 // Tmp array of jets form algoritm
761 Float_t etaAlgoJet[30];
762 Float_t phiAlgoJet[30];
763 Float_t etAlgoJet[30];
764 Int_t ncellsAlgoJet[30];
769 Int_t * index = new Int_t[nCell];
770 TMath::Sort(nCell, etCell, index);
772 // Variable used in centroide loop
790 for(Int_t icell = 0; icell < nCell; icell++)
792 Int_t jcell = index[icell];
793 if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
794 if(flagCell[jcell] != 0) continue; // if cell was used before
796 eta = etaCell[jcell];
797 phi = phiCell[jcell];
808 for(Int_t kcell =0; kcell < nCell; kcell++)
810 Int_t lcell = index[kcell];
811 if(lcell == jcell) continue; // cell itself
812 if(flagCell[lcell] != 0) continue; // cell used before
813 if(etCell[lcell] > etCell[jcell]) continue; // can this happen
815 deta = etaCell[lcell] - eta;
816 dphi = TMath::Abs(phiCell[lcell] - phi);
817 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
818 dr = TMath::Sqrt(deta * deta + dphi * dphi);
820 // calculate offset from initiate cell
821 deta = etaCell[lcell] - eta0;
822 dphi = phiCell[lcell] - phi0;
823 if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
824 if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
825 etas = etas + etCell[lcell]*deta;
826 phis = phis + etCell[lcell]*dphi;
827 ets = ets + etCell[lcell];
828 //new weighted eta and phi including this cell
829 eta = eta0 + etas/ets;
830 phi = phi0 + phis/ets;
831 // if cone does not move much, just go to next step
832 dphib = TMath::Abs(phi - phib);
833 if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
834 dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
835 if(dr <= minmove) break;
836 // cone should not move more than max_mov
837 dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
844 } else { // store this loop information
852 }//end of cells loop looking centroide
854 //avoid cones overloap (to be implemented in the future)
856 //flag cells in Rc, estimate total energy in cone
857 Float_t etCone = 0.0;
860 rc = header->GetRadius();
862 for(Int_t ncell =0; ncell < nCell; ncell++)
864 if(flagCell[ncell] != 0) continue; // cell used before
866 deta = etaCell[ncell] - eta;
867 // if(deta <= rc){ // Added to improve velocity -> to be tested
868 dphi = phiCell[ncell] - phi;
869 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
870 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
871 // if(dphi <= rc){ // Added to improve velocity -> to be tested
872 dr = TMath::Sqrt(deta * deta + dphi * dphi);
873 if(dr <= rc){ // cell in cone
874 flagCell[ncell] = -1;
875 etCone+=etCell[ncell];
879 // } // end deta <= rc
880 // } // end dphi <= rc
883 // select jets with et > background
884 // estimate max fluctuation of background in cone
885 Double_t ncellin = (Double_t)nCellIn;
886 Double_t ntcell = (Double_t)nCell;
887 Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/(ntcell));
889 Double_t etcmin = etCone ; // could be used etCone - etmin !!
890 //decisions !! etbmax < etcmin
892 for(Int_t mcell =0; mcell < nCell; mcell++)
894 if(flagCell[mcell] == -1){
896 flagCell[mcell] = 1; //flag cell as used
898 flagCell[mcell] = 0; // leave it free
901 //store tmp jet info !!!
904 etaAlgoJet[nJets] = eta;
905 phiAlgoJet[nJets] = phi;
906 etAlgoJet[nJets] = etCone;
907 ncellsAlgoJet[nJets] = nCellIn;
911 } // end of cells loop
913 for(Int_t p = 0; p < nJets; p++)
915 etaJet[p] = etaAlgoJet[p];
916 phiJet[p] = phiAlgoJet[p];
917 etJet[p] = etAlgoJet[p];
918 etallJet[p] = etAlgoJet[p];
919 ncellsJet[p] = ncellsAlgoJet[p];
927 ////////////////////////////////////////////////////////////////////////
928 void AliUA1JetFinderV2::RunAlgoritmC(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
929 Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
930 Float_t* etallJet, Int_t* ncellsJet)
933 // Check enough space! *to be done*
934 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
935 Float_t etCell[60000]; //! Cell Energy
936 Float_t etaCell[60000]; //! Cell eta
937 Float_t phiCell[60000]; //! Cell phi
938 Int_t flagCell[60000]; //! Cell flag
941 TAxis* xaxis = fLego->GetXaxis();
942 TAxis* yaxis = fLego->GetYaxis();
944 for (Int_t i = 1; i <= header->GetLegoNbinEta(); i++)
946 for (Int_t j = 1; j <= header->GetLegoNbinPhi(); j++)
948 e = fLego->GetBinContent(i,j);
949 if (e < 0.0) continue; // don't include this cells
950 Float_t eta = xaxis->GetBinCenter(i);
951 Float_t phi = yaxis->GetBinCenter(j);
953 etaCell[nCell] = eta;
954 phiCell[nCell] = phi;
955 flagCell[nCell] = 0; //default
960 // Parameters from header
961 Float_t minmove = header->GetMinMove();
962 Float_t maxmove = header->GetMaxMove();
963 Float_t rc = header->GetRadius();
964 Float_t etseed = header->GetEtSeed();
966 // Tmp array of jets form algoritm
967 Float_t etaAlgoJet[30];
968 Float_t phiAlgoJet[30];
969 Float_t etAlgoJet[30];
970 Int_t ncellsAlgoJet[30];
975 Int_t * index = new Int_t[nCell];
976 TMath::Sort(nCell, etCell, index);
977 // variable used in centroide loop
995 for(Int_t icell = 0; icell < nCell; icell++)
997 Int_t jcell = index[icell];
998 if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
999 if(flagCell[jcell] != 0) continue; // if cell was used before
1001 eta = etaCell[jcell];
1002 phi = phiCell[jcell];
1007 ets = etCell[jcell];
1013 for(Int_t kcell =0; kcell < nCell; kcell++)
1015 Int_t lcell = index[kcell];
1016 if(lcell == jcell) continue; // cell itself
1017 if(flagCell[lcell] != 0) continue; // cell used before
1018 if(etCell[lcell] > etCell[jcell]) continue; // can this happen
1020 deta = etaCell[lcell] - eta;
1021 dphi = TMath::Abs(phiCell[lcell] - phi);
1022 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1023 dr = TMath::Sqrt(deta * deta + dphi * dphi);
1026 // calculate offset from initiate cell
1027 deta = etaCell[lcell] - eta0;
1028 dphi = phiCell[lcell] - phi0;
1029 if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
1030 if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
1031 etas = etas + etCell[lcell]*deta;
1032 phis = phis + etCell[lcell]*dphi;
1033 ets = ets + etCell[lcell];
1034 //new weighted eta and phi including this cell
1035 eta = eta0 + etas/ets;
1036 phi = phi0 + phis/ets;
1037 // if cone does not move much, just go to next step
1038 dphib = TMath::Abs(phi - phib);
1039 if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
1040 dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
1041 if(dr <= minmove) break;
1042 // cone should not move more than max_mov
1043 dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
1050 } else { // store this loop information
1058 }//end of cells loop looking centroide
1060 // Avoid cones overloap (to be implemented in the future)
1062 // Flag cells in Rc, estimate total energy in cone
1063 Float_t etCone = 0.0;
1066 rc = header->GetRadius();
1067 for(Int_t ncell =0; ncell < nCell; ncell++)
1069 if(flagCell[ncell] != 0) continue; // cell used before
1071 deta = etaCell[ncell] - eta;
1072 dphi = phiCell[ncell] - phi;
1073 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1074 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1075 dr = TMath::Sqrt(deta * deta + dphi * dphi);
1076 if(dr <= rc){ // cell in cone
1077 flagCell[ncell] = -1;
1078 etCone+=etCell[ncell];
1084 // Select jets with et > background
1085 // estimate max fluctuation of background in cone
1086 Double_t ncellin = (Double_t)nCellIn;
1087 Double_t ntcell = (Double_t)nCell;
1088 Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
1090 Double_t etcmin = etCone ; // could be used etCone - etmin !!
1091 //decisions !! etbmax < etcmin
1093 for(Int_t mcell =0; mcell < nCell; mcell++){
1094 if(flagCell[mcell] == -1){
1096 flagCell[mcell] = 1; //flag cell as used
1098 flagCell[mcell] = 0; // leave it free
1101 //store tmp jet info !!!
1103 if(etbmax < etcmin) {
1104 etaAlgoJet[nJets] = eta;
1105 phiAlgoJet[nJets] = phi;
1106 etAlgoJet[nJets] = etCone;
1107 ncellsAlgoJet[nJets] = nCellIn;
1111 } // end of cells loop
1113 //reorder jets by et in cone
1114 //sort jets by energy
1115 Int_t * idx = new Int_t[nJets];
1116 TMath::Sort(nJets, etAlgoJet, idx);
1117 for(Int_t p = 0; p < nJets; p++)
1119 etaJet[p] = etaAlgoJet[idx[p]];
1120 phiJet[p] = phiAlgoJet[idx[p]];
1121 etJet[p] = etAlgoJet[idx[p]];
1122 etallJet[p] = etAlgoJet[idx[p]];
1123 ncellsJet[p] = ncellsAlgoJet[idx[p]];
1132 ////////////////////////////////////////////////////////////////////////
1133 void AliUA1JetFinderV2::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, const Float_t* ptT,
1134 const Int_t*vectT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1135 const Float_t* cFlag2T, const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet,
1136 const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet)
1139 // Background subtraction using cone method but without correction in dE/deta distribution
1140 // Cases to take into account the EMCal geometry are included
1143 //calculate energy inside and outside cones
1144 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1145 Int_t fOpt = fReader->GetReaderHeader()->GetDetector();
1146 fDebug = fReader->GetReaderHeader()->GetDebug();
1147 Float_t rc= header->GetRadius();
1151 for(Int_t j=0;j<30;j++){etIn[j]=0.;}
1153 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1155 for(Int_t ijet=0; ijet<nJ; ijet++){
1157 Float_t deta = etaT[jpart] - etaJet[ijet];
1158 Float_t dphi = phiT[jpart] - phiJet[ijet];
1159 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1160 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1162 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1163 if(dr <= rc){ // particles inside this cone
1164 multJet[ijet]+=vectT[jpart];
1165 injet[jpart] = ijet;
1167 if(cFlagT[jpart] == 1 || cFlag2T[jpart] == 1){ // pt cut
1168 etIn[ijet] += ptT[jpart];
1169 if(sFlagT[jpart] == 1) etsigJet[ijet]+= ptT[jpart];
1175 if(injet[jpart] == -1 && (cFlagT[jpart] == 1 || cFlag2T[jpart] == 1)){
1176 etOut += ptT[jpart]; // particle outside cones and pt cut
1178 } //end particle loop
1180 //estimate jets and background areas
1182 if(fOpt == 0 || fOpt == 1){
1183 Float_t areaJet[30];
1184 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1186 for(Int_t k=0; k<nJ; k++){
1187 Float_t detamax = etaJet[k] + rc;
1188 Float_t detamin = etaJet[k] - rc;
1189 Float_t accmax = 0.0; Float_t accmin = 0.0;
1190 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1191 Float_t h = header->GetLegoEtaMax() - etaJet[k];
1192 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1194 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1195 Float_t h = header->GetLegoEtaMax() + etaJet[k];
1196 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1198 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1199 areaOut = areaOut - areaJet[k];
1201 //subtract background using area method
1202 for(Int_t ljet=0; ljet<nJ; ljet++){
1203 Float_t areaRatio = areaJet[ljet]/areaOut;
1204 etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1207 // estimate new total background
1208 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1209 etbgTotalN = etOut*areaT/areaOut;
1211 else { // If EMCal included
1212 Float_t areaJet[30];
1213 Float_t areaOut = 2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax() - header->GetLegoPhiMin());
1214 for(Int_t k=0; k<nJ; k++){
1215 Float_t detamax = etaJet[k] + rc;
1216 Float_t detamin = etaJet[k] - rc;
1217 Float_t dphimax = phiJet[k] + rc;
1218 Float_t dphimin = phiJet[k] - rc;
1219 Float_t eMax = header->GetLegoEtaMax();
1220 Float_t eMin = header->GetLegoEtaMin();
1221 Float_t pMax = header->GetLegoPhiMax();
1222 Float_t pMin = header->GetLegoPhiMin();
1223 Float_t accetamax = 0.0; Float_t accetamin = 0.0;
1224 Float_t accphimax = 0.0; Float_t accphimin = 0.0;
1225 if((detamax > eMax && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1226 (detamax > eMax && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1227 Float_t h = eMax - etaJet[k];
1228 accetamax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1230 if((detamin < eMin && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1231 (detamin < eMin && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1232 Float_t h = eMax + etaJet[k];
1233 accetamin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1235 if((dphimax > pMax && detamax >= (eMin+2*rc) && detamax <= eMax )||
1236 (dphimax > pMax && detamin <= (eMax-2*rc) && detamin >= eMin )){
1237 Float_t h = pMax - phiJet[k];
1238 accphimax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1240 if((dphimin < eMin && detamax >= (eMin+2*rc) && detamax <= eMax )||
1241 (dphimin < eMin && detamin <= (eMax-2*rc) && detamin >= eMin )){
1242 Float_t h = phiJet[k] - pMin;
1243 accphimin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1246 if(detamax > eMax && dphimax > pMax ){
1247 Float_t he = eMax - etaJet[k];
1248 Float_t hp = pMax - phiJet[k];
1249 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1250 Float_t alphae = TMath::ACos(he/rc);
1251 Float_t alphap = TMath::ACos(hp/rc);
1252 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1254 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1255 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1258 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1259 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1260 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1261 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1265 if(detamax > eMax && dphimin < pMin ){
1266 Float_t he = eMax - etaJet[k];
1267 Float_t hp = phiJet[k] - pMin;
1268 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1269 Float_t alphae = TMath::ACos(he/rc);
1270 Float_t alphap = TMath::ACos(hp/rc);
1271 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1273 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1274 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1277 accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1278 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1279 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1280 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1284 if(detamin < eMin && dphimax > pMax ){
1285 Float_t he = eMax + etaJet[k];
1286 Float_t hp = pMax - phiJet[k];
1287 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1288 Float_t alphae = TMath::ACos(he/rc);
1289 Float_t alphap = TMath::ACos(hp/rc);
1290 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1292 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1293 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1296 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1297 accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1298 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1299 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1303 if(detamin < eMin && dphimin < pMin ){
1304 Float_t he = eMax + etaJet[k];
1305 Float_t hp = phiJet[k] - pMin;
1306 Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1307 Float_t alphae = TMath::ACos(he/rc);
1308 Float_t alphap = TMath::ACos(hp/rc);
1309 Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1311 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1312 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1315 accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1316 accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1317 ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1318 rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1321 areaJet[k] = rc*rc*TMath::Pi() - accetamax - accetamin - accphimax - accphimin;
1322 areaOut = areaOut - areaJet[k];
1323 } // end loop on jets
1325 //subtract background using area method
1326 for(Int_t ljet=0; ljet<nJ; ljet++){
1327 Float_t areaRatio = areaJet[ljet]/areaOut;
1328 etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1331 // estimate new total background
1332 Float_t areaT = 2*(header->GetLegoEtaMax()*header->GetLegoPhiMax());
1333 etbgTotalN = etOut*areaT/areaOut;
1338 ////////////////////////////////////////////////////////////////////////
1339 void AliUA1JetFinderV2::SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN,
1340 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
1341 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1342 Int_t* multJet, Int_t* injet)
1344 //background subtraction using cone method but without correction in dE/deta distribution
1346 //calculate energy inside and outside cones
1347 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1348 Float_t rc= header->GetRadius();
1351 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1352 // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
1353 for(Int_t ijet=0; ijet<nJ; ijet++){
1354 Float_t deta = etaT[jpart] - etaJet[ijet];
1355 Float_t dphi = phiT[jpart] - phiJet[ijet];
1356 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1357 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1358 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1359 if(dr <= rc){ // particles inside this cone
1361 injet[jpart] = ijet;
1362 if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
1363 etIn[ijet] += ptT[jpart];
1364 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
1369 if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
1370 etOut += ptT[jpart]; // particle outside cones and pt cut
1371 } //end particle loop
1373 //estimate jets and background areas
1374 Float_t areaJet[30];
1375 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1376 for(Int_t k=0; k<nJ; k++){
1377 Float_t detamax = etaJet[k] + rc;
1378 Float_t detamin = etaJet[k] - rc;
1379 Float_t accmax = 0.0; Float_t accmin = 0.0;
1380 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1381 Float_t h = header->GetLegoEtaMax() - etaJet[k];
1382 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1384 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1385 Float_t h = header->GetLegoEtaMax() + etaJet[k];
1386 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1388 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1389 areaOut = areaOut - areaJet[k];
1391 //subtract background using area method
1392 for(Int_t ljet=0; ljet<nJ; ljet++){
1393 Float_t areaRatio = areaJet[ljet]/areaOut;
1394 etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1397 // estimate new total background
1398 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1399 etbgTotalN = etOut*areaT/areaOut;
1404 ////////////////////////////////////////////////////////////////////////
1405 void AliUA1JetFinderV2::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
1406 const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1407 const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
1408 Float_t* etsigJet, Int_t* multJet, Int_t* injet)
1411 //background subtraction using statistical method
1412 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1413 Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
1415 //calculate energy inside
1416 Float_t rc= header->GetRadius();
1419 for(Int_t jpart = 0; jpart < nIn; jpart++)
1420 { // loop for all particles in array
1422 for(Int_t ijet=0; ijet<nJ; ijet++)
1424 Float_t deta = etaT[jpart] - etaJet[ijet];
1425 Float_t dphi = phiT[jpart] - phiJet[ijet];
1426 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1427 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1428 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1429 if(dr <= rc){ // particles inside this cone
1431 injet[jpart] = ijet;
1432 if(cFlagT[jpart] == 1){ // pt cut
1433 etIn[ijet]+= ptT[jpart];
1434 if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1439 } //end particle loop
1442 Float_t areaJet[30];
1443 Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1444 for(Int_t k=0; k<nJ; k++)
1446 Float_t detamax = etaJet[k] + rc;
1447 Float_t detamin = etaJet[k] - rc;
1448 Float_t accmax = 0.0; Float_t accmin = 0.0;
1449 if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1450 Float_t h = header->GetLegoEtaMax() - etaJet[k];
1451 accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1453 if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1454 Float_t h = header->GetLegoEtaMax() + etaJet[k];
1455 accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1457 areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1460 //subtract background using area method
1461 for(Int_t ljet=0; ljet<nJ; ljet++){
1462 Float_t areaRatio = areaJet[ljet]/areaOut;
1463 etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
1466 etbgTotalN = etbgStat;
1469 ////////////////////////////////////////////////////////////////////////
1470 void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
1471 Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1472 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1473 Int_t* multJet, Int_t* injet)
1475 // Cone background subtraction method taking into acount dEt/deta distribution
1476 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1478 Float_t rc= header->GetRadius();
1479 Float_t etamax = header->GetLegoEtaMax();
1480 Float_t etamin = header->GetLegoEtaMin();
1483 // jet energy and area arrays
1486 for(Int_t mjet=0; mjet<nJ; mjet++){
1487 char hEtname[256]; char hAreaname[256];
1488 sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1489 hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
1490 hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
1492 // background energy and area
1493 TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
1494 TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
1497 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1498 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1499 Float_t deta = etaT[jpart] - etaJet[ijet];
1500 Float_t dphi = phiT[jpart] - phiJet[ijet];
1501 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1502 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1503 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1504 if(dr <= rc){ // particles inside this cone
1505 injet[jpart] = ijet;
1507 if(cFlagT[jpart] == 1){// pt cut
1508 hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
1509 if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1515 if(injet[jpart] == -1 && cFlagT[jpart] == 1)
1516 hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
1517 } //end particle loop
1520 Float_t eta0 = etamin;
1521 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1522 Float_t eta1 = eta0 + etaw;
1523 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1524 Float_t etac = eta0 + etaw/2.0;
1525 Float_t areabg = etaw*2.0*TMath::Pi();
1526 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1527 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1528 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1529 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1530 Float_t areaj = 0.0;
1531 if(deta0 > rc && deta1 < rc){
1532 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1535 if(deta0 < rc && deta1 > rc){
1536 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1539 if(deta0 < rc && deta1 < rc){
1540 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1541 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1542 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
1543 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1544 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1546 hAreaJet[ijet]->Fill(etac,areaj);
1547 areabg = areabg - areaj;
1549 hAreaBackg->Fill(etac,areabg);
1552 } // end loop for all eta bins
1554 //subtract background
1555 for(Int_t kjet=0; kjet<nJ; kjet++){
1556 etJet[kjet] = 0.0; // first clear etJet for this jet
1557 for(Int_t bin = 0; bin< ndiv; bin++){
1558 if(hAreaJet[kjet]->GetBinContent(bin)){
1559 Float_t areab = hAreaBackg->GetBinContent(bin);
1560 Float_t etb = hEtBackg->GetBinContent(bin);
1561 Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1562 etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
1567 // calc background total
1568 Double_t etOut = hEtBackg->Integral();
1569 Double_t areaOut = hAreaBackg->Integral();
1570 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1571 etbgTotalN = etOut*areaT/areaOut;
1574 for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
1575 delete hEtJet[ljet];
1576 delete hAreaJet[ljet];
1583 ////////////////////////////////////////////////////////////////////////
1584 void AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
1585 Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1586 Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1587 Int_t* multJet, Int_t* injet)
1589 // Ratio background subtraction method taking into acount dEt/deta distribution
1590 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1591 //factor F calc before
1592 Float_t bgRatioCut = header->GetBackgCutRatio();
1595 Float_t rc= header->GetRadius();
1596 Float_t etamax = header->GetLegoEtaMax();
1597 Float_t etamin = header->GetLegoEtaMin();
1600 // jet energy and area arrays
1603 for(Int_t mjet=0; mjet<nJ; mjet++){
1604 char hEtname[256]; char hAreaname[256];
1605 sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1606 hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax); // change range
1607 hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax); // change range
1609 // background energy and area
1610 TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); // change range
1611 TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); // change range
1614 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1615 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1616 Float_t deta = etaT[jpart] - etaJet[ijet];
1617 Float_t dphi = phiT[jpart] - phiJet[ijet];
1618 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1619 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1620 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1621 if(dr <= rc){ // particles inside this cone
1623 injet[jpart] = ijet;
1624 if(cFlagT[jpart] == 1){ //pt cut
1625 hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
1626 if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1631 if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
1632 } //end particle loop
1635 Float_t eta0 = etamin;
1636 Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1637 Float_t eta1 = eta0 + etaw;
1638 for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1639 Float_t etac = eta0 + etaw/2.0;
1640 Float_t areabg = etaw*2.0*TMath::Pi();
1641 for(Int_t ijet=0; ijet<nJ; ijet++){ // loop for all jets
1642 Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1643 Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1644 Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1645 Float_t areaj = 0.0;
1646 if(deta0 > rc && deta1 < rc){
1647 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1650 if(deta0 < rc && deta1 > rc){
1651 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1654 if(deta0 < rc && deta1 < rc){
1655 acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1656 acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1657 if(eta1<etaJet[ijet]) areaj = acc1-acc0; // case 1
1658 if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1659 if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1661 hAreaJet[ijet]->Fill(etac,areaj);
1662 areabg = areabg - areaj;
1664 hAreaBackg->Fill(etac,areabg);
1667 } // end loop for all eta bins
1669 //subtract background
1670 for(Int_t kjet=0; kjet<nJ; kjet++){
1671 etJet[kjet] = 0.0; // first clear etJet for this jet
1672 for(Int_t bin = 0; bin< ndiv; bin++){
1673 if(hAreaJet[kjet]->GetBinContent(bin)){
1674 Float_t areab = hAreaBackg->GetBinContent(bin);
1675 Float_t etb = hEtBackg->GetBinContent(bin);
1676 Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1677 etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
1682 // calc background total
1683 Double_t etOut = hEtBackg->Integral();
1684 Double_t areaOut = hAreaBackg->Integral();
1685 Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1686 etbgTotalN = etOut*areaT/areaOut;
1689 for(Int_t ljet=0; ljet<nJ; ljet++){ // loop for all jets
1690 delete hEtJet[ljet];
1691 delete hAreaJet[ljet];
1698 ////////////////////////////////////////////////////////////////////////
1699 void AliUA1JetFinderV2::Reset()
1703 AliJetFinder::Reset();
1706 ////////////////////////////////////////////////////////////////////////
1707 void AliUA1JetFinderV2::WriteJHeaderToFile()
1709 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1713 ////////////////////////////////////////////////////////////////////////
1714 void AliUA1JetFinderV2::InitTask(TChain* tree)
1717 // initializes some variables
1718 AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1720 fLego = new TH2F("legoH","eta-phi",
1721 header->GetLegoNbinEta(), header->GetLegoEtaMin(),
1722 header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
1723 header->GetLegoPhiMin(), header->GetLegoPhiMax());
1725 fDebug = fReader->GetReaderHeader()->GetDebug();
1726 fOpt = fReader->GetReaderHeader()->GetDetector();
1728 // Tasks initialization
1730 fReader->CreateTasks(tree);