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