Adding include path to allow compilation of CleanGeom task
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALEMCGeometry.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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: AliEMCALEMCGeometry.cxx 29514 2008-10-26 10:24:38Z hristov $*/
17
18 //_________________________________________________________________________
19 // Geometry class  for EMCAL : singleton  
20 // EMCAL consists of layers of scintillator and lead
21 // with scintillator fiber arranged as "shish-kebab" skewers 
22 // Places the the Barrel Geometry of The EMCAL at Midrapidity
23 // between 80 and 180(or 190) degrees of Phi and
24 // -0.7 to 0.7 in eta 
25 //
26 //     EMCAL geometry tree:
27 //     EMCAL -> superModule -> module -> tower(cell)
28 //     Indexes
29 //     absId -> nSupMod     -> nModule -> (nIphi,nIeta)
30 //
31 //   Name choices: 
32 //   EMCAL_PDC06 (geometry used for PDC06 simulations, kept for backward compatibility)
33 //      = equivalent to SHISH_77_TRD1_2X2_FINAL_110DEG in old notation
34 //   EMCAL_COMPLETE (geometry for expected complete detector)
35 //      = equivalent to SHISH_77_TRD1_2X2_FINAL_110DEG scTh=0.176 pbTh=0.144
36 //          in old notation
37 //
38 //   EMCAL_FIRSTYEAR - geometry for December 2009 to December 2010 run period
39 //                     with four Super Modules
40 //  
41 //   Adding V1 (EMCAL_FIRSTYEARV1, EMCAL_COMPLETEV1) - geometry from December 2009 ; 
42 //                1. Fixed bug for positions of modules inside SM
43 //                   (first module has tilt 0.75 degree);
44 //                2. Added Al front plate (width 1 cm) and 2 paper sheets per sampling
45 //                   layer (additional 0.2 mm) 
46 //                   The sizes have updated with last information from production
47 //                   drawing (end of October 2010). 
48 //                3. COMPLETEV1 contains now only 10 SM for runs for year 2011
49 //                4. COMPLETE12SMV1 contains 12 SM for runs from year 2012 and on
50 //                5. COMPLETE12SMV1_DCAL contains 12 SM and 6 DCal SM
51 //                6. COMPLETE12SMV1_DCAL_DEV contains 12 SM and 10 DCal SM
52 //                7. COMPLETE12SMV1_DCAL_8SM contains 12 SM and 6 DCal SM and 2 extentions
53 //
54 //   EMCAL_WSUC (Wayne State test stand)
55 //      = no definite equivalent in old notation, was only used by
56 //          Aleksei, but kept for testing purposes
57 //
58 //   etc.
59 //
60 //
61 //
62 //*-- Author: Sahal Yacoob (LBL / UCT)
63 //     and  : Yves Schutz (SUBATECH)
64 //     and  : Jennifer Klay (LBL)
65 //     and  : Aleksei Pavlinov (WSU) 
66 //     and  : Magali Estienne (SUBATECH)
67 //     and  : Adapted for DCAL by M.L. Wang CCNU Wuhan & Subatech Oct-23-2009
68
69 // --- Root header files ---
70 #include <TObjArray.h>
71 #include <TObjString.h>
72 #include <TRegexp.h>
73
74 // -- ALICE Headers.
75 #include "AliLog.h"
76
77 // --- EMCAL headers
78 #include "AliEMCALEMCGeometry.h"
79 #include <cassert>
80
81 ClassImp(AliEMCALEMCGeometry)
82
83 // these initialisations are needed for a singleton
84 Bool_t    AliEMCALEMCGeometry::fgInit      = kFALSE;
85 const Char_t*   AliEMCALEMCGeometry::fgkDefaultGeometryName = "EMCAL_COMPLETE12SMV1";
86
87
88 AliEMCALEMCGeometry::AliEMCALEMCGeometry() 
89   : TNamed(),
90     fGeoName(0),fEMCSMSystem(0x0),fArrayOpts(0),fNAdditionalOpts(0),fECPbRadThickness(0.),fECScintThick(0.),
91     fNECLayers(0),fArm1PhiMin(0.),fArm1PhiMax(0.),fArm1EtaMin(0.),fArm1EtaMax(0.),fIPDistance(0.),
92     fShellThickness(0.),fZLength(0.),fDCALInnerEdge(0.),fDCALPhiMin(0),fDCALPhiMax(0),fEMCALPhiMax(0),
93     fDCALStandardPhiMax(0),fDCALInnerExtandedEta(0),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
94     fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
95     fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fPhiSuperModule(0),fNPhiSuperModule(0),
96     fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInModule(0),
97     // Trigger staff
98     fNTRUEta(0), fNTRUPhi(0), fNModulesInTRUEta(0), fNModulesInTRUPhi(0), fNEtaSubOfTRU(0), fNTotalTRU(0),
99     // 
100     fTrd1Angle(0.),f2Trd1Dx2(0.),fPhiGapForSM(0.),fKey110DEG(0),fnSupModInDCAL(0),fPhiBoundariesOfSM(0),
101     fPhiCentersOfSM(0),fPhiCentersOfSMSec(0),fEtaMaxOfTRD1(0),fTrd1AlFrontThick(0.0), fTrd1BondPaperThick(0.),
102     fCentersOfCellsEtaDir(0), fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),
103     fEtaCentersOfCells(0),fPhiCentersOfCells(0),fShishKebabTrd1Modules(0),
104     fParSM(), fILOSS(-1), fIHADR(-1),
105     //obsolete member data
106      fGap2Active(0.), fSteelFrontThick(0.), fTrd2AngleY(0.),
107     f2Trd2Dy2(0.), fEmptySpace(0.), fTubsR(0.), fTubsTurnAngle(0.)
108
109   // default ctor only for internal usage (singleton)
110   // must be kept public for root persistency purposes, 
111   // but should never be called by the outside world    
112   fParSM[0]=0; fParSM[1]=0; fParSM[2]=0;
113   fEnvelop[0] = 0; fEnvelop[1] = 0; fEnvelop[2] = 0;
114   for(Int_t i = 0; i < 6; i++) fkAdditionalOpts[i] = "";
115   
116   AliDebug(2, "AliEMCALEMCGeometry : default ctor ");
117 }
118 //______________________________________________________________________
119 AliEMCALEMCGeometry::AliEMCALEMCGeometry(const Text_t* name, const Text_t* title,
120                                          const Text_t* mcname, const Text_t* mctitle ) :
121   TNamed(name,title),
122     fGeoName(0),fEMCSMSystem(0x0),fArrayOpts(0),fNAdditionalOpts(0),fECPbRadThickness(0.),fECScintThick(0.),
123     fNECLayers(0),fArm1PhiMin(0.),fArm1PhiMax(0.),fArm1EtaMin(0.),fArm1EtaMax(0.),fIPDistance(0.),
124     fShellThickness(0.),fZLength(0.),fDCALInnerEdge(0.),fDCALPhiMin(0),fDCALPhiMax(0),fEMCALPhiMax(0),
125     fDCALStandardPhiMax(0),fDCALInnerExtandedEta(0),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
126     fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
127     fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fPhiSuperModule(0),fNPhiSuperModule(0),
128     fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInModule(0),
129     // Trigger staff
130     fNTRUEta(0), fNTRUPhi(0), fNModulesInTRUEta(0), fNModulesInTRUPhi(0), fNEtaSubOfTRU(0), fNTotalTRU(0),
131     // 
132     fTrd1Angle(0.),f2Trd1Dx2(0.),fPhiGapForSM(0.),fKey110DEG(0),fnSupModInDCAL(0),fPhiBoundariesOfSM(0),
133     fPhiCentersOfSM(0),fPhiCentersOfSMSec(0),fEtaMaxOfTRD1(0),fTrd1AlFrontThick(0.0), fTrd1BondPaperThick(0.),
134     fCentersOfCellsEtaDir(0),fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),
135     fEtaCentersOfCells(0),fPhiCentersOfCells(0),fShishKebabTrd1Modules(0),
136     fParSM(),fILOSS(-1), fIHADR(-1), 
137     //obsolete member data
138     fGap2Active(0.), fSteelFrontThick(0.), fTrd2AngleY(0.),
139     f2Trd2Dy2(0.), fEmptySpace(0.), fTubsR(0.), fTubsTurnAngle(0.)
140 {
141   // ctor only for internal usage (singleton)
142   AliDebug(2, Form("AliEMCALEMCGeometry(%s,%s,%s,%s) ", name,title,mcname,mctitle));
143
144   Init(mcname,mctitle);
145
146   //  CreateListOfTrd1Modules();
147
148   if (AliDebugLevel()>=2) {
149     PrintGeometry();
150   }
151
152 }
153 //______________________________________________________________________
154 AliEMCALEMCGeometry::AliEMCALEMCGeometry(const AliEMCALEMCGeometry& geom)
155   : TNamed(geom),
156     fGeoName(geom.fGeoName),
157     fEMCSMSystem(geom.fEMCSMSystem),
158     fArrayOpts(geom.fArrayOpts),
159     fNAdditionalOpts(geom.fNAdditionalOpts),
160     fECPbRadThickness(geom.fECPbRadThickness),
161     fECScintThick(geom.fECScintThick),
162     fNECLayers(geom.fNECLayers),
163     fArm1PhiMin(geom.fArm1PhiMin),
164     fArm1PhiMax(geom.fArm1PhiMax),
165     fArm1EtaMin(geom.fArm1EtaMin),
166     fArm1EtaMax(geom.fArm1EtaMax),
167     fIPDistance(geom.fIPDistance),
168     fShellThickness(geom.fShellThickness),
169     fZLength(geom.fZLength),
170     fDCALInnerEdge(geom.fDCALInnerEdge),
171     fDCALPhiMin(geom.fDCALPhiMin),
172     fDCALPhiMax(geom.fDCALPhiMax),
173     fEMCALPhiMax(geom.fEMCALPhiMax),
174     fDCALStandardPhiMax(geom.fDCALStandardPhiMax),
175     fDCALInnerExtandedEta(geom.fDCALInnerExtandedEta),
176     fNZ(geom.fNZ),
177     fNPhi(geom.fNPhi),
178     fSampling(geom.fSampling),
179     fNumberOfSuperModules(geom.fNumberOfSuperModules),
180     fFrontSteelStrip(geom.fFrontSteelStrip),
181     fLateralSteelStrip(geom.fLateralSteelStrip),
182     fPassiveScintThick(geom.fPassiveScintThick),
183     fPhiModuleSize(geom.fPhiModuleSize),
184     fEtaModuleSize(geom.fEtaModuleSize),
185     fPhiTileSize(geom.fPhiTileSize),
186     fEtaTileSize(geom.fEtaTileSize),
187     fLongModuleSize(geom.fLongModuleSize),
188     fPhiSuperModule(geom.fPhiSuperModule),
189     fNPhiSuperModule(geom.fNPhiSuperModule),
190     fNPHIdiv(geom.fNPHIdiv),
191     fNETAdiv(geom.fNETAdiv),
192     fNCells(geom.fNCells),
193     fNCellsInSupMod(geom.fNCellsInSupMod),
194     fNCellsInModule(geom.fNCellsInModule),
195     // Trigger staff
196     fNTRUEta(geom.fNTRUEta),
197     fNTRUPhi(geom.fNTRUPhi),
198     fNModulesInTRUEta(geom.fNModulesInTRUEta),
199     fNModulesInTRUPhi(geom.fNModulesInTRUPhi),
200     fNEtaSubOfTRU(geom.fNEtaSubOfTRU),
201     fNTotalTRU(geom.fNTotalTRU),
202     //
203     fTrd1Angle(geom.fTrd1Angle),
204     f2Trd1Dx2(geom.f2Trd1Dx2),
205     fPhiGapForSM(geom.fPhiGapForSM),
206     fKey110DEG(geom.fKey110DEG),
207     fnSupModInDCAL(geom.fnSupModInDCAL),
208     fPhiBoundariesOfSM(geom.fPhiBoundariesOfSM),
209     fPhiCentersOfSM(geom.fPhiCentersOfSM),
210     fPhiCentersOfSMSec(geom.fPhiCentersOfSMSec),
211     fEtaMaxOfTRD1(geom.fEtaMaxOfTRD1),
212     fTrd1AlFrontThick(geom.fTrd1AlFrontThick),
213     fTrd1BondPaperThick(geom.fTrd1BondPaperThick),
214     fCentersOfCellsEtaDir(geom.fCentersOfCellsEtaDir),
215     fCentersOfCellsXDir(geom.fCentersOfCellsXDir),
216     fCentersOfCellsPhiDir(geom.fCentersOfCellsPhiDir),
217     fEtaCentersOfCells(geom.fEtaCentersOfCells),
218     fPhiCentersOfCells(geom.fPhiCentersOfCells),
219     fShishKebabTrd1Modules(geom.fShishKebabTrd1Modules),
220     fILOSS(geom.fILOSS), fIHADR(geom.fIHADR),
221     //obsolete member data
222     fGap2Active(geom.fGap2Active),
223     fSteelFrontThick(geom.fSteelFrontThick),
224     fTrd2AngleY(geom.fTrd2AngleY),
225     f2Trd2Dy2(geom.f2Trd2Dy2),
226     fEmptySpace(geom.fEmptySpace),
227     fTubsR(geom.fTubsR),
228     fTubsTurnAngle(geom.fTubsTurnAngle)
229 {
230   //copy ctor
231   fParSM[0]=geom.fParSM[0]; 
232   fParSM[1]=geom.fParSM[1]; 
233   fParSM[2]=geom.fParSM[2];
234   fEnvelop[0] = geom.fEnvelop[0]; 
235   fEnvelop[1] = geom.fEnvelop[1]; 
236   fEnvelop[2] = geom.fEnvelop[2];
237   for(Int_t i = 0; i < 6; i++) fkAdditionalOpts[i] = geom.fkAdditionalOpts[i];
238
239 }
240
241 //______________________________________________________________________
242 AliEMCALEMCGeometry::~AliEMCALEMCGeometry(void){
243     // dtor
244  delete fEMCSMSystem;
245 }
246
247 //______________________________________________________________________
248 void AliEMCALEMCGeometry::Init(const Text_t* mcname, const Text_t* mctitle){
249   //
250   // Initializes the EMCAL parameters based on the name
251   // Only Shashlyk geometry is available, but various combinations of
252   // layers and number of supermodules can be selected with additional
253   // options or geometry name
254   //
255
256   fkAdditionalOpts[0] = "nl=";       // number of sampling layers (fNECLayers)
257   fkAdditionalOpts[1] = "pbTh=";     // cm, Thickness of the Pb   (fECPbRadThick)
258   fkAdditionalOpts[2] = "scTh=";     // cm, Thickness of the Sc    (fECScintThick)
259   fkAdditionalOpts[3] = "latSS=";    // cm, Thickness of lateral steel strip (fLateralSteelStrip)
260   fkAdditionalOpts[4] = "allILOSS="; // = 0,1,2,3,4 (4 - energy loss without fluctuation)
261   fkAdditionalOpts[5] = "allIHADR="; // = 0,1,2 (0 - no hadronic interaction)
262
263   fNAdditionalOpts = sizeof(fkAdditionalOpts) / sizeof(char*);
264
265   // geometry
266   fgInit = kFALSE; // Assume failed until proven otherwise.
267   fGeoName   = GetName();
268   fGeoName.ToUpper();
269
270   //Convert old geometry names to new ones
271   if(fGeoName.Contains("SHISH_77_TRD1_2X2_FINAL_110DEG")) {
272     if(fGeoName.Contains("PBTH=0.144") && fGeoName.Contains("SCTH=0.176")) {
273       fGeoName = "EMCAL_COMPLETE";
274     } else {
275       fGeoName = "EMCAL_PDC06";
276     }
277   }
278   
279   if(fGeoName.Contains("WSUC")) fGeoName = "EMCAL_WSUC";
280
281   //check that we have a valid geometry name
282   if(!( fGeoName.Contains("EMCAL_PDC06")
283      || fGeoName.Contains("EMCAL_WSUC")  
284      || fGeoName.Contains("EMCAL_COMPLETE")
285      || fGeoName.Contains("EMCAL_COMPLETEV1")
286      || fGeoName.Contains("EMCAL_COMPLETE12SMV1") 
287      || fGeoName.Contains("EMCAL_FIRSTYEAR")
288      || fGeoName.Contains("EMCAL_FIRSTYEARV1") )) {
289     Fatal("Init", "%s is an undefined geometry!", fGeoName.Data()) ; 
290   }
291
292   // Option to know whether we have the "half" supermodule(s) or not
293   fKey110DEG = 0;
294   if(  fGeoName.Contains("COMPLETE")
295     || fGeoName.Contains("PDC06") 
296     || fGeoName.Contains("12SM")) fKey110DEG = 1; // for GetAbsCellId
297   if(fGeoName.Contains("COMPLETEV1"))  fKey110DEG = 0; 
298   fShishKebabTrd1Modules = 0;
299
300   fnSupModInDCAL = 0;
301   if(fGeoName.Contains("DCAL_DEV")){
302     fnSupModInDCAL = 10;  
303   } else if(fGeoName.Contains("DCAL_8SM")){
304     fnSupModInDCAL = 8;
305   } else if(fGeoName.Contains("DCAL")){
306     fnSupModInDCAL = 6;
307   }
308
309   // JLK 13-Apr-2008
310   //default parameters are those of EMCAL_COMPLETE geometry
311   //all others render variations from these at the end of
312   //geometry-name specific options
313
314   fNumberOfSuperModules = 12;       // 12 = 6 * 2 (6 in phi, 2 in Z)
315   fNPhi                 = 12;       // module granularity in phi within smod (azimuth)
316   fNZ                   = 24;       // module granularity along Z within smod (eta)
317   fNPHIdiv = fNETAdiv   = 2;        // tower granularity within module
318   fArm1PhiMin           = 80.0;     // degrees, Starting EMCAL Phi position
319   fArm1PhiMax           = 200.0;    // degrees, Ending EMCAL Phi position
320   fArm1EtaMin           = -0.7;     // pseudorapidity, Starting EMCAL Eta position
321   fArm1EtaMax           = +0.7;     // pseudorapidity, Ending EMCAL Eta position
322   fIPDistance           = 428.0;    // cm, radial distance to front face from nominal vertex point
323   fPhiGapForSM          = 2.;       // cm, only for final TRD1 geometry
324   fFrontSteelStrip      = 0.025;    // 0.025cm = 0.25mm  (13-may-05 from V.Petrov)
325   fPassiveScintThick    = 0.8;      // 0.8cm   = 8mm     (13-may-05 from V.Petrov)
326   fLateralSteelStrip    = 0.01;     // 0.01cm  = 0.1mm   (13-may-05 from V.Petrov) - was 0.025
327   fTrd1Angle            = 1.5;      // in degrees       
328
329   fSampling             = 1.;       // should be calculated with call to DefineSamplingFraction()
330   fNECLayers            = 77;       // (13-may-05 from V.Petrov) - can be changed with additional options
331   fECScintThick         = 0.176;    // scintillator layer thickness
332   fECPbRadThickness     = 0.144;    // lead layer thickness
333
334   fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
335   fEtaModuleSize = fPhiModuleSize;
336
337   fZLength              = 700.;     // Z coverage (cm)
338   fPhiSuperModule       = 20. ;     // phi in degree
339   fDCALInnerEdge        = fIPDistance * TMath::Tan( fTrd1Angle * 8.* TMath::DegToRad());
340
341   //needs to be called for each geometry and before setting geometry
342   //parameters which can depend on the outcome
343   CheckAdditionalOptions();
344
345   //modifications to the above for PDC06 geometry
346   if(fGeoName.Contains("PDC06")){ // 18-may-05 - about common structure
347     fECScintThick  = fECPbRadThickness = 0.16;// (13-may-05 from V.Petrov)    
348     CheckAdditionalOptions();
349   }
350
351   //modifications to the above for WSUC geometry
352   if(fGeoName.Contains("WSUC")){ // 18-may-05 - about common structure
353     fNumberOfSuperModules = 2;  // 27-may-05; Nov 24,2010 for TB
354     fNPhi = fNZ = 4; 
355     fTrd1AlFrontThick   = 1.0;  // one cm
356     // Bond paper - two sheets around Sc tile
357     fTrd1BondPaperThick = 0.01; // 0.01cm = 0.1 mm
358     
359     fPhiModuleSize = 12.0;
360     fEtaModuleSize = fPhiModuleSize;
361     fLateralSteelStrip = 0.015; // 0.015cm  = 0.15mm
362
363     CheckAdditionalOptions();
364   }
365
366   //In 2009-2010 data taking runs only 4 SM, in the upper position.
367   if(fGeoName.Contains("FIRSTYEAR")){   
368     fNumberOfSuperModules = 4;  
369     fArm1PhiMax           = 120.0;
370     CheckAdditionalOptions();   
371   }
372   
373   if(fGeoName.Contains("FIRSTYEARV1") || fGeoName.Contains("COMPLETEV1") || fGeoName.Contains("COMPLETE12SMV1") ){
374     // Oct 26,2010 : First module has tilt = 0.75 degree : 
375     // look to AliEMCALShishKebabTrd1Module::DefineFirstModule(key)
376     // New sizes from production drawing, added Al front plate.
377     // The thickness of sampling is change due to existing two sheets of paper.
378     
379     // Will replace fFrontSteelStrip
380     fTrd1AlFrontThick   = 1.0;  // one cm
381     // Bond paper - two sheets around Sc tile
382     fTrd1BondPaperThick = 0.01; // 0.01cm = 0.1 mm
383     
384     fPhiModuleSize = 12.0;
385     fEtaModuleSize = fPhiModuleSize;
386     fLateralSteelStrip = 0.015; // 0.015cm  = 0.15mm
387     
388     if(fGeoName.Contains("COMPLETEV1"))
389     {
390       fNumberOfSuperModules = 10;       
391       fArm1PhiMax           = 180.0;
392     }
393     else if (fGeoName.Contains("COMPLETE12SMV1"))
394     {
395       fNumberOfSuperModules = 12;       
396       fArm1PhiMax           = 200.0; 
397     }
398     if (fGeoName.Contains("DCAL"))
399     {
400       fNumberOfSuperModules = 12 + fnSupModInDCAL;
401       fArm1PhiMax           = 320.0; 
402       if(fGeoName.Contains("DCAL_8SM"))      fArm1PhiMax      = 340.0;      // degrees, End of DCAL Phi position
403       else if(fGeoName.Contains("DCAL_DEV")) fArm1PhiMin      = 40.0;       // degrees, Starting EMCAL(shifted) Phi position
404       fDCALPhiMin           = fArm1PhiMax - 10.*fnSupModInDCAL; 
405     }
406     CheckAdditionalOptions();   
407   }
408
409    fEMCSMSystem = new Int_t[fNumberOfSuperModules];
410    Int_t iSM = 0;
411
412  // BASIC EMCAL SM
413    if(fGeoName.Contains("WSUC") ){
414      for(int i = 0; i<2; i++){
415        fEMCSMSystem[iSM] = kEMCAL_Standard;
416        iSM++;
417      }
418    } else if(fGeoName.Contains("FIRSTYEAR") ){
419      for(int i = 0; i<4; i++){
420        fEMCSMSystem[iSM] = kEMCAL_Standard;
421        iSM++;
422      }
423    } else if( fGeoName.Contains("PDC06")
424            || fGeoName.Contains("COMPLETE") ){
425      for(int i = 0; i<10; i++){
426        fEMCSMSystem[iSM] = kEMCAL_Standard;
427        iSM++;
428      }
429    }
430  // EMCAL 110SM
431    if(fKey110DEG && fGeoName.Contains("12SM") ){
432      for(int i = 0; i<2; i++){
433        fEMCSMSystem[iSM] = kEMCAL_Half;
434        if(fGeoName.Contains("12SMV1") ){
435          fEMCSMSystem[iSM] = kEMCAL_3rd;
436        }
437        iSM++;
438      }
439    }
440  // DCAL SM
441    if(fnSupModInDCAL && fGeoName.Contains("DCAL")){
442      if(fGeoName.Contains("8SM")) {
443        for(int i = 0; i<fnSupModInDCAL-2; i++){
444          fEMCSMSystem[iSM] = kDCAL_Standard;
445          iSM++;
446        }
447        for(int i = 0; i<2; i++){
448          fEMCSMSystem[iSM] = kDCAL_Ext;
449          iSM++;
450        }
451      } else {
452        for(int i = 0; i<fnSupModInDCAL; i++){
453          fEMCSMSystem[iSM] = kDCAL_Standard;
454          iSM++;
455        }
456      }
457    }
458
459   // constant for transition absid <--> indexes
460   fNCellsInModule = fNPHIdiv*fNETAdiv;
461   fNCellsInSupMod = fNCellsInModule*fNPhi*fNZ;
462   fNCells = 0;
463    for( int i=0; i<fNumberOfSuperModules; i++) {
464      if(      GetSMType(i) == kEMCAL_Standard) fNCells +=   fNCellsInSupMod   ;
465      else if( GetSMType(i) == kEMCAL_Half)     fNCells +=   fNCellsInSupMod/2 ;
466      else if( GetSMType(i) == kEMCAL_3rd)      fNCells +=   fNCellsInSupMod/3 ;
467      else if( GetSMType(i) == kDCAL_Standard)  fNCells += 2*fNCellsInSupMod/3 ;
468      else if( GetSMType(i) == kDCAL_Ext)       fNCells +=   fNCellsInSupMod/3 ;
469      else AliError(Form("Uknown SuperModule Type !!"));
470    }
471
472   fNPhiSuperModule = fNumberOfSuperModules/2;
473   if(fNPhiSuperModule < 1) fNPhiSuperModule = 1;
474     
475   fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05 
476   fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05 
477
478   fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);  
479   if(fGeoName.Contains("V1")){
480     Double_t ws = fECScintThick + fECPbRadThickness + 2.*fTrd1BondPaperThick; // sampling width
481     // Number of Pb tiles = Number of Sc tiles - 1
482     fLongModuleSize = fTrd1AlFrontThick + (ws*fNECLayers - fECPbRadThickness);
483   }
484   f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
485
486   if(!fGeoName.Contains("WSUC")) fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
487
488   //These parameters are used to create the mother volume to hold the supermodules
489   //2cm padding added to allow for misalignments - JLK 30-May-2008
490   fEnvelop[0]     = fIPDistance - 1.; // mother volume inner radius
491   fEnvelop[1]     = fIPDistance + fShellThickness + 1.; // mother volume outer r.
492   fEnvelop[2]     = fZLength + 2.; //mother volume length 
493
494   // Local coordinates
495   fParSM[0] = GetShellThickness()/2.;        
496   fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
497   fParSM[2] = fZLength/4.;  //divide by 4 to get half-length of SM
498
499   // SM phi boundaries - (0,1),(2,3) ... - has the same boundaries;
500   fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
501   fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
502   fPhiCentersOfSMSec.Set(fNumberOfSuperModules/2);
503   Double_t kfSupermodulePhiWidth = fPhiSuperModule*TMath::DegToRad();
504   fPhiCentersOfSM[0]    = (fArm1PhiMin + fPhiSuperModule/2.) * TMath::DegToRad(); // Define from First SM
505   fPhiCentersOfSMSec[0] = fPhiCentersOfSM[0];  // the same in the First SM
506   fPhiBoundariesOfSM[0] = fPhiCentersOfSM[0] - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
507   fPhiBoundariesOfSM[1] = fPhiCentersOfSM[0] + TMath::ATan2(fParSM[1] , fIPDistance);
508   if(fNumberOfSuperModules > 2) { // 2 to Max
509     Int_t tmpSMType = GetSMType(2);
510      for(int i = 1; i<fNPhiSuperModule; i++) {
511        fPhiBoundariesOfSM[2*i]   += fPhiBoundariesOfSM[2*i-2] + kfSupermodulePhiWidth;
512        if(tmpSMType == GetSMType(2*i)) {
513          fPhiBoundariesOfSM[2*i+1]  += fPhiBoundariesOfSM[2*i-1] + kfSupermodulePhiWidth;
514        } else { 
515          //changed SM Type, redefine the [2*i+1] Boundaries
516          tmpSMType = GetSMType(2*i);
517          if(        GetSMType(2*i)  == kEMCAL_Standard) {
518            fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[2*i] + kfSupermodulePhiWidth;
519          } else if( GetSMType(2*i)  == kEMCAL_Half)     {
520            fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[2*i] + 2.*TMath::ATan2((fParSM[1])/2, fIPDistance);
521          } else if( GetSMType(2*i)  == kEMCAL_3rd )     {
522            fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[2*i] + 2.*TMath::ATan2((fParSM[1])/3, fIPDistance);
523          } else if( GetSMType(2*i)  == kDCAL_Standard ) {      // jump the gap
524            fPhiBoundariesOfSM[2*i]   = (fDCALPhiMin - fArm1PhiMin)*TMath::DegToRad() + fPhiBoundariesOfSM[0];
525            fPhiBoundariesOfSM[2*i+1] = (fDCALPhiMin - fArm1PhiMin)*TMath::DegToRad() + fPhiBoundariesOfSM[1];
526          } else if( GetSMType(2*i)  == kDCAL_Ext)       {
527            fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[2*i] + 2.*TMath::ATan2((fParSM[1])/3, fIPDistance); 
528          }
529       }
530       fPhiCentersOfSM[i]    = (fPhiBoundariesOfSM[2*i]+fPhiBoundariesOfSM[2*i+1])/2.;
531       fPhiCentersOfSMSec[i] = fPhiBoundariesOfSM[2*i] + TMath::ATan2(fParSM[1] , fIPDistance);
532     }
533   }
534
535   //inner extend in eta (same as outer part) for DCal (0.189917), //calculated from the smallest gap (1# cell to the 80-degree-edge),
536   Double_t fInnerExtandedPhi = 1.102840997; //calculated from the smallest gap (1# cell to the 80-degree-edge), too complicatd to explain...
537   fDCALInnerExtandedEta = -TMath::Log(TMath::Tan( (TMath::Pi()/2. - 8*fTrd1Angle*TMath::DegToRad() + (TMath::Pi()/2 - fNZ*fTrd1Angle*TMath::DegToRad() - TMath::ATan(TMath::Exp(fArm1EtaMin))*2))/2.));
538
539   fEMCALPhiMax  = fArm1PhiMin;    
540   fDCALPhiMax   = fDCALPhiMin;// DCAl extention will not be included
541   for( Int_t i = 0; i < fNumberOfSuperModules; i+=2) {
542     if(      GetSMType(i) == kEMCAL_Standard ) fEMCALPhiMax += 20.;
543     else if( GetSMType(i) == kEMCAL_Half     ) fEMCALPhiMax += fPhiSuperModule/2. + fInnerExtandedPhi;
544     else if( GetSMType(i) == kEMCAL_3rd      ) fEMCALPhiMax += fPhiSuperModule/3. + 4.0*fInnerExtandedPhi/3.0;
545     else if( GetSMType(i) == kDCAL_Standard  ) {fDCALPhiMax  += 20.; fDCALStandardPhiMax = fDCALPhiMax;}
546     else if( GetSMType(i) == kDCAL_Ext       ) fDCALPhiMax  += fPhiSuperModule/3. + 4.0*fInnerExtandedPhi/3.0;
547     else AliError("Unkown SM Type!!");
548   }
549 // for compatible reason
550 // if(fNumberOfSuperModules == 4) {fEMCALPhiMax = fArm1PhiMax ;}
551  if(fNumberOfSuperModules == 12) {fEMCALPhiMax = fArm1PhiMax ;}  
552   
553   //called after setting of scintillator and lead layer parameters
554   DefineSamplingFraction(mcname,mctitle);
555
556   
557   // TRU parameters - Apr 29,08 by PAI. 
558   // These parameters values was updated at Nov 05, 2007
559   // As is on Olivier  BOURRION (LPSC) ppt preasentation 
560   // at ALICE trigger meeting at 13th-14th March
561   fNTRUEta = 1;           // was 3
562   fNTRUPhi = 3;           // was 1
563   fNModulesInTRUEta = 24; // was 8
564   fNModulesInTRUPhi = 4;  // was 12
565   // Jet trigger 
566   // 3*6*10 + 2*6*2 = 204 -> matrix (nphi(17), neta(12))
567   fNEtaSubOfTRU     = 6;  
568
569   fNTotalTRU = 0;
570   for(Int_t i = 0; i < GetNumberOfSuperModules(); i++){
571     if(      GetSMType(i) == kEMCAL_Standard)  fNTotalTRU += 3;
572     else if( GetSMType(i) == kEMCAL_Half)      fNTotalTRU += 1;
573     else if( GetSMType(i) == kEMCAL_3rd)       fNTotalTRU += 1;
574     else if( GetSMType(i) == kDCAL_Standard)   fNTotalTRU += 3;
575     else if( GetSMType(i) == kDCAL_Ext)        fNTotalTRU += 1;
576     else {
577       AliError(Form("Uknown SuperModule Type !!"));
578     }
579   }
580
581   fgInit = kTRUE; 
582 }
583
584 //___________________________________________________________________
585 void AliEMCALEMCGeometry::PrintGeometry()
586 {
587   // Separate routine is callable from broswer; Nov 7,2006
588   printf("\nInit: geometry of EMCAL named %s :\n", fGeoName.Data());
589   if(fArrayOpts) {
590     for(Int_t i=0; i<fArrayOpts->GetEntries(); i++){
591       TObjString *o = (TObjString*)fArrayOpts->At(i);
592       printf(" %i : %s \n", i, o->String().Data());
593     }
594   }
595   if(fGeoName.Contains("DCAL")) {
596    printf("Phi min of DCAL SuperModule: %7.1f, DCAL has %d SuperModule\n", fDCALPhiMin, fnSupModInDCAL);
597    printf("The DCAL inner edge is +- %7.1f\n", fDCALInnerEdge);
598     if(fGeoName.Contains("DCAL_8SM")) printf("DCAL has its 2 EXTENTION SM\n");
599   }
600   printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
601   printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f -> for EMCAL envelope only\n",  
602            GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
603
604   printf( "               ECAL      : %d x (%f cm Pb, %f cm Sc) \n", 
605   GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; 
606   printf("                fSampling %5.2f \n",  fSampling );
607   printf(" fIPDistance       %6.3f cm \n", fIPDistance);
608   printf(" fNPhi %i   |  fNZ %i \n", fNPhi, fNZ);
609   printf(" fNCellsInModule %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInModule, fNCellsInSupMod, fNCells);
610   printf(" X:Y module size     %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
611   printf(" X:Y   tile size     %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
612   printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
613   printf(" fLongModuleSize     %6.3f cm \n", fLongModuleSize);
614   printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
615   printf(" supermodule width in phi direction %f \n", fPhiSuperModule );
616   printf(" fILOSS %i : fIHADR %i \n", fILOSS, fIHADR);
617   printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
618   printf(" f2Trd1Dx2  %7.4f\n",  f2Trd1Dx2);
619   printf(" fTrd1AlFrontThick   %7.4f \n", fTrd1AlFrontThick);
620   printf(" fTrd1BondPaperThick %5.4f \n", fTrd1BondPaperThick);
621   printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", 
622          fParSM[0],fParSM[1],fParSM[2]);
623   printf(" fPhiGapForSM  %7.4f cm (%7.4f <- phi size in degree)\n",  
624          fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg());
625   if( fKey110DEG && !fGeoName.Contains("12SMV1") ) printf(" Last two modules have size 10 degree in  phi (180<phi<190)\n");
626   if( fKey110DEG && fGeoName.Contains("12SMV1")) printf(" Last two modules have size 6.6 degree in  phi (180<phi<186.6)\n");
627   printf(" phi SM boundaries \n"); 
628   for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
629     printf(" %i : %7.15f(%7.12f) -> %7.15f(%7.12f) : center %7.15f(%7.12f) \n", i, 
630            fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(),
631            fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(),
632            fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg());
633   }
634
635 }
636
637 //______________________________________________________________________
638 void AliEMCALEMCGeometry::CheckAdditionalOptions()
639 {
640   // Feb 06,2006
641   // Additional options that
642   // can be used to select
643   // the specific geometry of 
644   // EMCAL to run
645   // Dec 27,2006
646   // adeed allILOSS= and allIHADR= for MIP investigation
647   fArrayOpts = new TObjArray;
648   Int_t nopt = ParseString(fGeoName, *fArrayOpts);
649   if(nopt==1) { // no aditional option(s)
650     fArrayOpts->Delete();
651     delete fArrayOpts;
652     fArrayOpts = 0; 
653     return;
654   }              
655   for(Int_t i=1; i<nopt; i++){
656     TObjString *o = (TObjString*)fArrayOpts->At(i); 
657
658     TString addOpt = o->String();
659     Int_t indj=-1;
660     for(Int_t j=0; j<fNAdditionalOpts; j++) {
661       TString opt = fkAdditionalOpts[j];
662       if(addOpt.Contains(opt,TString::kIgnoreCase)) {
663           indj = j;
664         break;
665       }
666     }
667     if(indj<0) {
668       AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n", 
669                       addOpt.Data()));
670       assert(0);
671     } else {
672       AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n", 
673                       addOpt.Data(), indj, fkAdditionalOpts[indj]));
674       if       (addOpt.Contains("NL=",TString::kIgnoreCase))   {// number of sampling layers
675         sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
676         AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
677       } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
678         sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
679       } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
680         sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
681       } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
682         sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
683         AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
684       } else if(addOpt.Contains("ILOSS=",TString::kIgnoreCase)) {// As in Geant
685         sscanf(addOpt.Data(),"ALLILOSS=%i", &fILOSS);
686         AliDebug(2,Form(" fILOSS %i \n", fILOSS));
687       } else if(addOpt.Contains("IHADR=",TString::kIgnoreCase)) {// As in Geant
688         sscanf(addOpt.Data(),"ALLIHADR=%i", &fIHADR);
689         AliDebug(2,Form(" fIHADR %i \n", fIHADR));
690       }
691     }
692   }
693 }
694
695 //__________________________________________________________________
696 void AliEMCALEMCGeometry::DefineSamplingFraction(const Text_t* mcname, const Text_t* mctitle)
697 {
698   // Jun 05,2006
699   // Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
700   // Keep for compatibilty
701   //
702   
703   // Sampling factor for G3
704   fSampling = 10.87; // Default value - Nov 25,2010
705   if(fNECLayers == 69) {        // 10% layer reduction
706     fSampling = 12.55;
707   } else if(fNECLayers == 61) { // 20% layer reduction
708     fSampling = 12.80;
709   } else if(fNECLayers == 77) {
710     if(fGeoName.Contains("V1")){
711       fSampling = 10.87; //Adding paper sheets and cover plate; Nov 25,2010
712     } else if   (fECScintThick>0.159 && fECScintThick<0.161) { // original sampling fraction, equal layers
713       fSampling = 12.327; // fECScintThick = fECPbRadThickness = 0.160;
714     } else if (fECScintThick>0.175 && fECScintThick<0.177) { // 10% Pb thicknes reduction
715       fSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;
716     } else if(fECScintThick>0.191 && fECScintThick<0.193) { // 20% Pb thicknes reduction
717       fSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;
718     }
719   }
720
721   // Default sampling factor for G3, modify it for other transport model
722   TString mcName  = mcname;
723   TString mcTitle = mctitle;
724
725   Float_t samplingFactorTranportModel = 1. ;
726   if     (mcName.Contains("Geant3")) samplingFactorTranportModel = 1.;//0.988 // Do nothing
727   else if(mcName.Contains("Fluka") ) samplingFactorTranportModel = 1.; // To be set
728   else if(mcName.Contains("Geant4")){
729     if(mcTitle.Contains("EMV"))      samplingFactorTranportModel = 1.096; // 0.906, 0.896 (OPT)
730     else                             samplingFactorTranportModel = 0.86; // 1.15 (CHIPS), 1.149 (BERT), 1.147 (BERT_CHIPS) 
731   }      
732   
733   AliDebug(2,Form("MC modeler <%s>, Title <%s>: Sampling %f, model fraction with respect to G3 %f, final sampling %f \n",
734                mcName.Data(),mcTitle.Data(),fSampling,samplingFactorTranportModel,fSampling*samplingFactorTranportModel));
735
736   
737   fSampling*=samplingFactorTranportModel;
738   
739 }
740
741 //________________________________________________________________________________________________
742 Double_t AliEMCALEMCGeometry::GetPhiCenterOfSMSec(Int_t nsupmod) const
743 {
744   //returns center of supermodule in phi
745   int i = nsupmod/2;
746   return fPhiCentersOfSMSec[i];
747
748 }
749
750 //________________________________________________________________________________________________
751 Double_t AliEMCALEMCGeometry::GetPhiCenterOfSM(Int_t nsupmod) const
752 {
753   //returns center of supermodule in phi
754   int i = nsupmod/2;
755   return fPhiCentersOfSM[i];
756
757 }
758
759 //________________________________________________________________________________________________
760 Bool_t AliEMCALEMCGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
761 {
762   // 0<= nSupMod <=17; phi in rad
763   static int i;
764   if(nSupMod<0 || nSupMod >12+fnSupModInDCAL-1) return kFALSE;
765   i = nSupMod/2;
766   phiMin = (Double_t)fPhiBoundariesOfSM[2*i];
767   phiMax = (Double_t)fPhiBoundariesOfSM[2*i+1];
768   return kTRUE;
769 }
770
771 //________________________________________________________________________________________________
772 Bool_t AliEMCALEMCGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
773 {
774   // 0<= nPhiSec <=max; phi in rad
775   // 0;  gap boundaries between  0th&2th  | 1th&3th SM
776   // 1;  gap boundaries between  2th&4th  | 3th&5th SM
777   // 2;  gap boundaries between  4th&6th  | 5th&7th SM
778   // 3;  gap boundaries between  6th&8th  | 7th&9th SM
779   // 4;  gap boundaries between  8th&10th | 9th&11th SM
780   // 5;  gap boundaries between 10th&12th | 11h&13th SM
781   //             ...
782   if(nPhiSec<0 || nPhiSec >5+fnSupModInDCAL/2-1) return kFALSE;
783   phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
784   phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
785   return kTRUE;
786 }
787
788 //________________________________________________________________________________________________
789 int AliEMCALEMCGeometry::ParseString(const TString &topt, TObjArray &Opt)
790
791         //Parse string, does what? GCB 08/09
792         Ssiz_t begin, index, end, end2;
793         begin = index = end = end2 = 0;
794         TRegexp separator("[^ ;,\\t\\s/]+");
795         while ( (begin < topt.Length()) && (index != kNPOS) ) {
796                 // loop over given options
797                 index = topt.Index(separator,&end,begin);
798                 if (index >= 0 && end >= 1) {
799                         TString substring(topt(index,end));
800                         Opt.Add(new TObjString(substring.Data()));
801                 }
802                 begin += end+1;
803         }
804         return Opt.GetEntries();
805 }
806