]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliJetGrid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix of detail in the background subtraction method. (Rafael Diaz Valdes)
[u/mrichter/AliRoot.git] / JETAN / AliJetGrid.cxx
CommitLineData
4a01bb2c 1/**************************************************************************
2 * Copyright(c) 2001-2002, 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: AliJetGrid.cxx,v 1.0 07/09/2006 */
17
18//=========================================================================
19// *** 07 September 2006
20// This class allows to fill a really general grid in (eta,phi)
21// Two types of grid can be setted :
22// if fGrid==0 -> the complete acceptance of a rectangular detector acceptance is filled
23// if fGrid==1 -> the previous grid - an other rectangular detector acceptance is filled
24// A (eta,phi(rad)) position can be extracted from an index value without looping over all entries
25// An index position can be extracted from a (eta,phi(rad)) position without looping over all entries
26//
27// How to run it :
28// > AliJetGrid *grid = new AliJetGrid((NbinPhi-1),(NbinEta-1),PhiMin,PhiMax,EtaMin,EtaMax);
29// > grid->SetGridType(...); // 0 or 1 for the moment
30// > grid->SetMatrixIndexes();
31// > grid->SetIndexIJ();
32//
33// Author : magali.estienne@ires.in2p3.fr
34//=========================================================================
35
36// Standard headers
37#include <Riostream.h>
38// Root headers
39#include <TMatrixD.h>
40#include <TArrayD.h>
41#include <TArrayI.h>
42// AliRoot headers
43#include "AliJetGrid.h"
44
87f40625 45ClassImp(AliJetGrid)
4a01bb2c 46
4a01bb2c 47
b92e2ccf 48//__________________________________________________________
49AliJetGrid::AliJetGrid():
50 fGrid(0),
51 fNphi(0),
52 fNeta(0),
53 fPhi(0),
54 fEta(0),
55 fIndex(0),
56 fIndexI(0),
57 fIndexJ(0),
58 fPhiMin(0),
59 fPhiMax(0),
60 fEtaMin(0),
61 fEtaMax(0),
62 fEtaBinInTPCAcc(0),
63 fPhiBinInTPCAcc(0),
64 fEtaBinInEMCalAcc(0),
65 fPhiBinInEMCalAcc(0),
66 fNbinEta(0),
67 fNbinPhi(0),
68 fMaxPhi(0),
69 fMinPhi(0),
70 fMaxEta(0),
71 fMinEta(0),
72 fDebug(1)
73 {
4a01bb2c 74 // Default constructor
4a01bb2c 75}
76
77//__________________________________________________________
b92e2ccf 78AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax):
79 fGrid(0),
80 fNphi(nphi),
81 fNeta(neta),
82 fPhi(0),
83 fEta(0),
84 fIndex(0),
85 fIndexI(0),
86 fIndexJ(0),
87 fPhiMin(0),
88 fPhiMax(0),
89 fEtaMin(0),
90 fEtaMax(0),
91 fEtaBinInTPCAcc(0),
92 fPhiBinInTPCAcc(0),
93 fEtaBinInEMCalAcc(0),
94 fPhiBinInEMCalAcc(0),
95 fNbinEta(0),
96 fNbinPhi(0),
97 fMaxPhi(phiMax),
98 fMinPhi(phiMin),
99 fMaxEta(etaMax),
100 fMinEta(etaMin),
101 fDebug(1)
102 {
4a01bb2c 103 // Standard constructor
b92e2ccf 104 fPhi = new TArrayD(fNphi+1);
105 fEta = new TArrayD(fNeta+1);
4a01bb2c 106 fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1);
107 fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1);
b92e2ccf 108
4a01bb2c 109 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = (phiMax-phiMin)/fNphi*i+phiMin;
110 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = (etaMax-etaMin)/fNeta*i+etaMin;
b92e2ccf 111
4a01bb2c 112 if(fDebug > 3){
b92e2ccf 113 for(Int_t i=0; i<(fNphi+1); i++) cout << (*fPhi)[i] << endl;
114 for(Int_t i=0; i<(fNeta+1); i++) cout << (*fEta)[i] << endl;
4a01bb2c 115 }
b92e2ccf 116
4a01bb2c 117 fIndex = new TMatrixD(fNphi+1,fNeta+1);
b92e2ccf 118
4a01bb2c 119}
120
121//__________________________________________________________
122AliJetGrid::AliJetGrid(const AliJetGrid& grid):TNamed(grid) {
123
124 // Copy constructor
125
126 fNphi = grid.fNphi;
127 fNeta = grid.fNeta;
128 fPhiMin = grid.fPhiMin;
129 fPhiMax = grid.fPhiMax;
130 fEtaMin = grid.fEtaMin;
131 fEtaMax = grid.fEtaMax;
132 fEtaBinInTPCAcc = grid.fEtaBinInTPCAcc;
133 fPhiBinInTPCAcc = grid.fPhiBinInTPCAcc;
134 fEtaBinInEMCalAcc = grid.fEtaBinInEMCalAcc;
135 fPhiBinInEMCalAcc = grid.fPhiBinInEMCalAcc;
136 fNbinPhi = grid.fNbinPhi;
137 fMaxPhi = grid.fMaxPhi;
138 fMinPhi = grid.fMinPhi;
139 fMaxEta = grid.fMaxEta;
140 fMinEta = grid.fMinEta;
141
142 fPhi = new TArrayD(fNphi+1);
143 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = grid.fPhi->At(i);
144 fEta = new TArrayD(fNeta+1);
145 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = grid.fEta->At(i);
146
147 fIndex = new TMatrixD(fNphi+1,fNeta+1);
148 for(Int_t i=0; i<fNphi+1; i++) {
149 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*grid.fIndex)(i,j);
150 }
151}
152
153//__________________________________________________________
154AliJetGrid::~AliJetGrid() {
155
156 // Destructor
157 delete fPhi;
158 delete fEta;
159 delete fIndexI;
160 delete fIndexJ;
161 delete fIndex;
162}
163
164//__________________________________________________________
165void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal)
166{ // To set initial parameters
167
168 fPhiMin = phiMinCal; // rad
169 fPhiMax = phiMaxCal; // rad
170 fEtaMin = etaMinCal;
171 fEtaMax = etaMaxCal;
172 fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
173
174 // Define some binning numbers
175 if(fGrid==0){
176 for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
177 fPhiBinInEMCalAcc = 0;
178
179 for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
180 fEtaBinInEMCalAcc = 0;
181 }
182
183 if(fGrid==1){
184 for(Int_t i=0; i<fNphi+1; i++) {
185 fPhiBinInTPCAcc++;
186 if(fPhi->At(i) >= fPhiMin &&
187 fPhi->At(i) <= fPhiMax)
188 fPhiBinInEMCalAcc++;
189 }
190 for(Int_t i=0; i<fNeta+1; i++) {
191 fEtaBinInTPCAcc++;
192 if((fEta->At(i) >= fEtaMin) &&
193 (fEta->At(i) <= fEtaMax))
194 fEtaBinInEMCalAcc++;
195 }
196 }
197
198}
199
200//__________________________________________________________
201TArrayD* AliJetGrid::GetArrayEta()
202{ // Returns an array with the eta points
203
204 return fEta;
205
206}
207
208//__________________________________________________________
209TArrayD* AliJetGrid::GetArrayPhi()
210{ // Returns an array with the phi points
211
212 return fPhi;
213
214}
215
216//__________________________________________________________
217TMatrixD* AliJetGrid::GetIndexObject()
218{ // Returns a pointer to the matrix
219
220 return fIndex;
221
222}
223
224//__________________________________________________________
225void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
226 Float_t &mineta, Float_t &maxeta)
227{ // Returns EMCAL acceptance initially setted
228
229 nphi = fNphi;
230 neta = fNeta;
231 minphi = fPhiMin;
232 maxphi = fPhiMax;
233 mineta = fEtaMin;
234 maxeta = fEtaMax;
235}
236
237//__________________________________________________________
238void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
239 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi)
240{ // Returns number of bins in TPC and EMCAL geometry
241
242 etabintpc = fEtaBinInTPCAcc;
243 phibintpc = fPhiBinInTPCAcc;
244 etabinemc = fEtaBinInEMCalAcc;
245 phibinemc = fPhiBinInEMCalAcc;
246 nbinphi = fNbinPhi;
247}
248
249//__________________________________________________________
250Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
251{ // Tells the index value of a corresponding (eta,phi) real position
252 // Loop over all entries -> takes time.
253 // Used one time at the begining to fill the grids
254
255 /* this is how bins are numbered
256
257 in all TPC or in TPC - EMCAL
258 ... ... ... ..
259 ---------------
260 ... ... . 10 | | | 11
261 ---+---+--- ---------------
262 ^ 6 | 7 | 8 or 8 | | | 9
263 | ---+---+--- ---------------
264 3 | 4 | 5 4 | 5 | 6 | 7
265 phi ---+---+--- ---------------
266 0 | 1 | 2 0 | 1 | 2 | 3
267
268
269 eta ->
270 */
271
272 Int_t etaBin=0,phiBin=0,absID=0;
273 Int_t etaBin2=0,etaBin3=0;
274
275 // Fill all the grid in eta/phi (all TPC acceptance)
276 //-----------------------------------------------------
277 if(fGrid==0){
278 if(eta <= fEta->At(0)) {
279 etaBin = 0;
280 } else if(eta >= fEta->At(fNeta)) {
281 etaBin = fNeta;
282 } else {
283 for(Int_t i=0; i<fNeta+1; i++) {
284 if(eta < fEta->At(i)) {
285 etaBin = i-1;
286 break;
287 }
288 }
289 }
290 if(phi <= fPhi->At(0)) {
291 phiBin = 0;
292 } else if(phi >= fPhi->At(fNphi)) {
293 phiBin = fNphi;
294 } else {
295 for(Int_t i=0; i<fNphi+1; i++) {
296 if(phi < fPhi->At(i)) {
297 phiBin = i-1;
298 break;
299 }
300 }
301 }
302
303 // Calculate absolute id
304 absID = phiBin*(fNeta+1) + etaBin;
305
306 }
307
308 // Fill the grid but do not count id in EMCal acceptance
309 //------------------------------------------------------
310 if((eta >= fEtaMin && eta <= fEtaMax) &&
311 (phi >= fPhiMin && phi <= fPhiMax)){
312 etaBin = etaBin2 = etaBin3 = -1;
313 phiBin = -1;
314 }
315
316 if(fGrid == 1){
317 if(phi<fPhiMin) {
318 if(eta <= fEta->At(0)) {
319 etaBin = 0;
320 } else if(eta >= fEta->At(fNeta)) {
321 etaBin = fNeta;
322 } else {
323 for(Int_t i=0; i<fNeta+1; i++) {
324 if(eta < fEta->At(i)) {
325 etaBin = i-1;
326 break;
327 }
328 }
329 }
330 }
331 else {
332 if((phi>=fPhiMin) && (phi<=fPhiMax)) {
333 if(eta <= fEta->At(0)) {
334 etaBin2 = 0;
335 } else if(eta >= fEta->At(fNeta)) {
336 etaBin2 = fNeta-fEtaBinInEMCalAcc;
337 } else {
338 for(Int_t i=0; i<fNeta+1; i++) {
339 if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) ||
340 (fEtaMax<eta && eta<fEta->At(fNeta)))) {
341 // cout << "i : " << i << endl;
342 if(eta<fEtaMin)
343 etaBin2 = i-1;
344 else etaBin2 = i-1-fEtaBinInEMCalAcc;
345 break;
346 }
347 }
348 }
349 }
350 else {
351 if(eta <= fEta->At(0)) {
352 etaBin3 = 0;
353 } else if(eta >= fEta->At(fNeta)) {
354 etaBin3 = fNeta;
355 } else {
356 for(Int_t i=0; i<fNeta+1; i++) {
357 if(eta < fEta->At(i)) {
358 etaBin3 = i-1;
359 break;
360 }
361 }
362 }
363 }
364 }
365
366 if(phi <= fPhi->At(0)) {
367 phiBin = 0;
368 } else if(phi >= fPhi->At(fNphi)) {
369 phiBin = fNphi;
370 } else {
371 for(Int_t i=0; i<fNphi+1; i++) {
372 if(phi < fPhi->At(i)) {
373 phiBin = i-1;
374 break;
375 }
376 }
377 }
378
379 // Calculate absID
380 //-------------------------------------------------------
381 if(phi<fPhiMin)
382 absID = phiBin*(fNeta+1) + etaBin;
383 if(phi>=fPhiMin && phi<=fPhiMax) {
384 if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
385 else{
386 absID = (fNbinPhi+1)*(fNeta+1)
387 + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
388 + etaBin2;
389 }
390 }
391 if(phi>fPhiMax)
392 absID = (fNbinPhi+1)*(fNeta+1)
393 + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
394 + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
395 + etaBin3;
396
397 } // END OPTION==1
398
399 return absID;
400
401}
402
403//__________________________________________________________
404void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
405{ // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
406
407 for(Int_t j=0; j<fNphi+1; j++) {
408 for(Int_t i=0; i<fNeta+1; i++) {
409
410 // TPC grid only
411 //-------------------------------------
412 if(fGrid==0) {
413 if(j*(fNeta+1)+i == index) {
414 eta = fEta->At(i);
415 phi = fPhi->At(j);
416 }
417 }
418
419 // TPC-EMCAL grid
420 //-------------------------------------
421 Int_t ii = 0;
422 if(i==0) ii = 0;
423 if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
424 if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
425
426 if(fGrid==1) {
427 if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
428 eta = fEta->At(i);
429 phi = fPhi->At(j);
430 }
431
432 if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
433 ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
434 if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
435 else eta = fEta->At(i);
436 phi = fPhi->At(j);
437 }
438
439 if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) &&
440 ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
441 +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
442 eta = fEta->At(i);
443 phi = fPhi->At(j);
444 }
445 }
446 }
447 }
448}
449
450//__________________________________________________________
451Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
452{ // Get index value for a (eta,phi) position - Direct value
453
454 Int_t ieta = GetIndexJFromEta(eta);
455 Int_t iphi = GetIndexIFromPhi(phi);
456
457 Int_t index = GetMatrixIndex(iphi,ieta);
458
459 if(fDebug>10){
460 cout << "(phi,eta) : " << phi << ", " << eta << endl;
461 cout << "index : " << index << endl;
462 }
463 return index;
464
465}
466
467//__________________________________________________________
468Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
469{ // Get eta id
470 // Eta discretized
471
472 Int_t idEta =0;
473 Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
474 if(fDebug>20)
475 {
476 cout << "eta : " << eta << endl;
477 cout << "fMaxEta : " << fMaxEta << endl;
478 cout << "fMinEta : " << fMinEta << endl;
479 cout << "fNeta : " << fNeta << endl;
480 cout << "temp eta before cast : " << temp << endl;
481 }
482 idEta = static_cast<Int_t>(temp+0.5);
483 if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
484 return idEta;
485}
486//__________________________________________________________
487Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
488{ // Get phi id
489 // Phi discretized
490
491 Int_t idPhi = 0;
492 Double_t temp = 0.;
493 if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
494 else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
495
496 if(fDebug>20)
497 {
498 cout << "phi : " << phi << endl;
499 cout << "fMaxPhi : " << fMaxPhi << endl;
500 cout << "fMinPhi : " << fMinPhi << endl;
501 cout << "fNphi : " << fNphi << endl;
502 cout << "temp phi before cast : " << temp << endl;
503 }
504 idPhi = static_cast<Int_t>(temp+0.5);
505 if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
506 return idPhi;
507
508
509}
510
511//__________________________________________________________
512void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
513{ // Allows to set parameters using only one index (if fGrid==0) !!
514 // Not used !
515
516 Int_t iphi = (Int_t)i/fNeta;
517 Int_t ieta = i-iphi*fNeta;
518 SetMatrixIndex(iphi,ieta,par);
519
520 return;
521}
522
523//__________________________________________________________
524void AliJetGrid::SetMatrixIndexes()
525{ // Fill the final matrix object with the corresponding index in eta/phi
526
527 for(Int_t i=0; i<fNphi+1; i++){
528 for(Int_t j=0; j<fNeta+1; j++){
529 (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
530 if(fDebug>2){
531 cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) <<
532 ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
533 }
534 }
535 }
536 printf("##########################################################\n");
537 printf("TMatrix object filled !\n");
538 printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
539 printf("##########################################################\n");
540}
541
542//__________________________________________________________
543void AliJetGrid::SetIndexIJ()
544{ //
545
546 for(Int_t i=0; i<fNphi+1; i++){
547 for(Int_t j=0; j<fNeta+1; j++){
548 Int_t id = static_cast<Int_t>((*fIndex)(i,j));
549
550 if(id!=-1)
551 {
552 (*fIndexI)[id] = i;
553 (*fIndexJ)[id] = j;
554 }
555 }
556 }
557
558 printf("##########################################################\n");
559 printf(" In SetIndexIJ - Grid indexes setted !\n");
560 printf("##########################################################\n");
561}
562
563//__________________________________________________________
564void AliJetGrid::GetIJFromIndex(Int_t index, Int_t i, Int_t j)
565{ // returns i position id of eta and j position id of phi for a given grid index
566 i = (*fIndexI)[index];
567 j = (*fIndexJ)[index];
568}
569
570//__________________________________________________________
571void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
572{ // returns eta, phi values for a given grid index
573
574 phi = fPhi->At((*fIndexI)[index]);
575 eta = fEta->At((*fIndexJ)[index]);
576}
577
578//__________________________________________________________
579Int_t AliJetGrid::GetNEntries()
580{ // Returns the number of entries of the grid
581
582 if(fDebug>20){
583 cout << "fMaxPhi : " << fMaxPhi << endl;
584 cout << "fMaxEta : " << fMaxEta << endl;
585 }
586
587 Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
588 if(fDebug>20) cout << "indexNum : " << indexNum << endl;
589 return indexNum;
590
591}
592
593//__________________________________________________________
594Int_t AliJetGrid::GetNEntries2()
595{ // Returns the number of entries of the grid
596
597 Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
598 if(fDebug>20) cout << "indexNum : " << indexNum << endl;
599 return indexNum;
600
601}
602
603
604
605
606
607