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