Federico: AliHLTEMCALGeometry handles transformation from local to global coordinates
[u/mrichter/AliRoot.git] / HLT / EMCAL / AliHLTEMCALGeometry.cxx
1 /**************************************************************************
2  * This file is property of and copyright by the Experimental Nuclear     *
3  * Physics Group, Dep. of Physics                                         *
4  * University of Oslo, Norway, 2007                                       *
5  *                                                                        *
6  * Author: federico ronchetti         for the ALICE HLT Project.*
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 #include "AliHLTEMCALGeometry.h"
18 #include "AliHLTEMCALConstants.h"
19 #include "AliHLTEMCALConstant.h"
20 #include "assert.h"
21 #include "AliHLTCaloConstantsHandler.h"
22 #include "AliHLTEMCALSharedMemoryInterface.h" 
23 #include "TVector3.h"
24
25 using namespace EmcalHLTConst;
26
27 ClassImp(AliHLTEMCALGeometry);
28 TGeoManager *gGeoManager = 0;
29
30 AliHLTEMCALGeometry::AliHLTEMCALGeometry() :
31         AliHLTCaloGeometry ("EMCAL"),
32         fShmPtr(0),
33         fGeo(0),
34         fEMCALGeometry(0)
35 {
36         //fGeo = new AliEMCALGeoUtils("EMCAL_COMPLETE","EMCAL");
37         //fGeo = new AliEMCALGeometry("EMCAL_COMPLETE","EMCAL");
38         //fGeo =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
39         //TGeoManager::Import("/home/fedro/work/AliRoot/test/QA/geometry.root");
40         fGeo = new AliEMCALGeoUtils("EMCAL_COMPLETE","EMCAL");
41         fShmPtr = new AliHLTEMCALSharedMemoryInterface();
42         GetGeometryFromCDB();
43 }
44
45
46 AliHLTEMCALGeometry::~AliHLTEMCALGeometry()
47 {
48
49 }
50   
51 void 
52 AliHLTEMCALGeometry::GetGlobalCoordinates(AliHLTEMCALRecPointDataStruct &recPoint, AliHLTCaloGlobalCoordinate &globalCoord)
53 {
54
55          Int_t istrip = 0;
56          Float_t z0 = 0;
57          Float_t zb = 0;
58          Float_t z_is = 0;
59          Float_t d = 0;
60
61          Float_t x,y,z; // return variables in terry's RF
62
63          Float_t dz = fCaloConstants->GetMINCELLSTEPETA(); // base cell width in eta
64          Float_t dx = fCaloConstants->GetCELLSTEPPHI(); // base cell width in phi
65
66           // parameters for shower depth calculation
67          Float_t X0 = fCaloConstants->GetRADLENGTH();
68          Float_t Ecr = fCaloConstants->GetCRITICENERGY();
69          Float_t Cj;
70
71          Float_t teta0 = fCaloConstants->GetCELLANGLE(); //tapering angle (deg)
72          Float_t teta1; //working angle
73
74          Float_t L = fCaloConstants->GetCELLHEIGHT();
75
76          // converting to MEV
77          Float_t E = recPoint.fAmp * 1000;
78
79          //TVector3 v1;
80                 Double_t glob[3];
81                 Double_t loc[3];
82
83          if (recPoint.fZ >= 47.5 || recPoint.fZ<-0.5) {
84                  Logging(kHLTLogError, "HLT", "EMCAL", "AliHLTEMCALGeometry::GetGlobalCoordinates: invalid Z recpoint ");
85                  return;
86           }
87
88           if (recPoint.fX >= 23.5 || recPoint.fX<-0.5) {
89                   Logging(kHLTLogError, "HLT", "EMCAL", "AliHLTEMCALGeometry::GetGlobalCoordinates: invalid X recpoint ");
90                   return;
91           }
92
93
94           switch ( recPoint.fParticle )
95           {
96           case 0:
97                   Cj = - fCaloConstants->GetCJ(); // photon
98                   d = X0 * TMath::Log( E / Ecr + Cj);
99                   break;
100
101           case 1:
102                   Cj = + fCaloConstants->GetCJ(); // electron
103                   d = X0 * TMath::Log( E / Ecr + Cj);
104                   break;
105
106           case 2:
107                   // hadron
108                   d = 0.5 * L;
109                   break;
110
111           default:
112                   Cj = - fCaloConstants->GetCJ(); // defaulting to photon
113                   d = X0 * TMath::Log( E / Ecr + Cj);
114           }
115
116            istrip = int ((recPoint.fZ + 0.5 ) / 2);
117
118            // module angle
119            teta1 = TMath::DegToRad() * istrip * teta0;
120
121            // calculation of module corner along z
122            // as a function of strip
123
124            for (Int_t is=0; is<= istrip; is++) {
125
126              teta1 = TMath::DegToRad() * is * teta0;
127
128              z_is = z_is + 2*dz*(TMath::Sin(teta1)*TMath::Tan(teta1) + TMath::Cos(teta1));
129
130            }
131
132            z0 = dz * (recPoint.fZ - 2*istrip + 0.5);
133            zb = (2*dz-z0-d*TMath::Tan(teta1))*TMath::Cos(teta1);
134
135            z = z_is - zb*TMath::Cos(teta1);
136
137 //         cout << "----> istrip: " << istrip << endl;
138 //         cout << "----> z0: "<< z0 << endl;
139 //         cout << "----> zb: "<< zb << endl;
140 //         cout << "----> corner z: "<< z_is << endl;
141
142 //         cout << "----> teta1: "<< TMath::RadToDeg()*teta1 << endl;
143
144            y = d/TMath::Cos(teta1) + zb*TMath::Sin(teta1);
145
146            x = (recPoint.fX + 0.5)*dx;
147
148            // cout << "x: " << x << " y: "<< y << " z " << z << endl;
149
150         // check coordinate origin
151         loc[0] = x;
152         loc[1] = y;
153         loc[2] = z;
154
155          if(!fGeo)
156    {
157       Logging(kHLTLogError, "HLT", "EMCAL", "AliHLTEMCALGeometry::GetGlobalCoordinates: no geometry initialised");
158       return;
159    }
160
161          ConvertRecPointCoordinates(loc[1], loc[2], loc[3]);
162
163          fGeo->GetGlobal(loc, glob, recPoint.fModule);
164
165          globalCoord.fX = glob[0];
166          globalCoord.fY = glob[1];
167          globalCoord.fZ = glob[2];
168 }
169  
170 void 
171 AliHLTEMCALGeometry::GetCellAbsId(UInt_t module, UInt_t y, UInt_t z, Int_t& AbsId)
172 {
173
174         if(!fGeo)
175         {
176                  Logging(kHLTLogError, "HLT", "EMCAL", "AliHLTEMCALGeometry::GetCellAbsId: no geometry initialised");
177                  return;
178
179         }
180         
181         AbsId = fGeo->GetAbsCellIdFromCellIndexes(module, y, z);
182         
183 }
184
185 void AliHLTEMCALGeometry::ConvertRecPointCoordinates(Double_t &x, Double_t &y, Double_t &z) const
186 {
187         Double_t DX = 13.869008;
188         Double_t DY = 72.559998;
189         Double_t DZ = 175.00000;
190
191          x = y - DX;  //fixme
192          y = -x + DY; //fixme
193          z = z - DZ;  //fixme
194
195 }
196
197
198 int
199 AliHLTEMCALGeometry::GetGeometryFromCDB()
200 {
201   cout << "Getting geometry..." << endl;
202
203  // AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
204
205   AliCDBPath path("GRP","Geometry","Data");
206   if(path.GetPath())
207     {
208       //      HLTInfo("configure from entry %s", path.GetPath());
209       AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path/*,GetRunNo()*/);
210       if (pEntry)
211         {
212           if(!fEMCALGeometry)
213             {
214               delete fEMCALGeometry;
215               fEMCALGeometry = 0;
216             }
217
218           gGeoManager = (TGeoManager*) pEntry->GetObject();
219           cout<< "gGeoManager = 0x%x" << gGeoManager << endl;
220
221           if(gGeoManager)
222             {
223               fGeo = new AliEMCALGeoUtils("EMCAL_COMPLETE","EMCAL");
224               //fClusterAnalyserPtr->SetGeometry(fEMCALGeometry);
225             }
226
227         }
228       else
229         {
230                   //HLTError("can not fetch object \"%s\" from OCDB", path);
231         }
232     }
233   return 0;
234 }