in last commit, had accidentally removed implementation of method PbInTrd1 which...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALv0.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 /* $Id$ */
16
17 //_________________________________________________________________________
18 // Implementation version v0 of EMCAL Manager class 
19 // An object of this class does not produce hits nor digits
20 // It is the one to use if you do not want to produce outputs in TREEH or TREED
21 // This class places a Geometry of the EMCAL in the ALICE Detector as defined in AliEMCALGeometry.cxx                 
22 //*-- Author: Yves Schutz (SUBATECH)
23 //*-- and   : Sahal Yacoob (LBL / UCT)
24 //          : Aleksei Pavlinov (WSU)     SHASHLYK
25
26 // --- ROOT system ---
27 #include <cassert>
28
29 #include <TNode.h>
30 #include <TBRIK.h>
31 #include <TTRD1.h>
32 #include <TTRAP.h>
33 #include <TPGON.h>
34 #include <TTUBS.h>
35 #include <TGeometry.h>
36 #include <TGeoPhysicalNode.h>
37 #include <TGeoManager.h>
38 #include <TGeoMatrix.h>
39 #include <TVirtualMC.h>
40 #include <TArrayI.h>
41 #include <TROOT.h>
42 #include <TList.h>
43 #include <TVector2.h>
44 #include <cassert>
45
46 #include "AliEMCALShishKebabModule.h"
47 #include "AliEMCALShishKebabTrd1Module.h"
48
49 // --- Standard library ---
50
51 //#include <stdio.h>
52
53 // --- AliRoot header files ---
54
55 #include "AliEMCALv0.h"
56 #include "AliEMCALGeometry.h"
57 #include "AliRun.h"
58 #include "AliLog.h"
59
60 ClassImp(AliEMCALv0)
61
62 // EMCAL material
63 enum
64  {
65   kIdAIR   = 1599, 
66   kIdPB    = 1600, 
67   kIdSC    = 1601, 
68   kIdSTEEL = 1603
69  };
70
71
72 //______________________________________________________________________
73 AliEMCALv0::AliEMCALv0()
74   : AliEMCAL(),
75     fShishKebabModules(),fEnvelop1(0),fIdRotm(0),fIdTmedArr(0),
76     fSampleWidth(0),fSmodPar0(0),fSmodPar1(0),fSmodPar2(0)
77 {
78   //default ctor
79 }
80
81 //______________________________________________________________________
82 AliEMCALv0::AliEMCALv0(const char *name, const char *title)
83   : AliEMCAL(name,title),
84     fShishKebabModules(),fEnvelop1(0),fIdRotm(0),fIdTmedArr(0),
85     fSampleWidth(0),fSmodPar0(0),fSmodPar1(0),fSmodPar2(0)
86 {
87   // ctor : title is used to identify the layout
88   // Apr 25, 2006
89   // Nov 22, 2006 - case of 1X1  
90   TString ntmp(GetTitle());
91   ntmp.ToUpper();
92   if(ntmp == "TRD1") { // TRD1 is alias for SHISH_77_TRD1_2X2_FINAL_110DEG NL=69
93     ntmp = "SHISH_77_TRD1_2X2_FINAL_110DEG NL=69";
94     SetTitle(ntmp.Data());
95   }
96
97   AliEMCALGeometry *geom = GetGeometry() ; 
98   //geom->CreateListOfTrd1Modules(); 
99   fShishKebabModules = geom->GetShishKebabTrd1Modules(); 
100 }
101
102 //______________________________________________________________________
103 void AliEMCALv0::BuildGeometry()
104 {
105     // Display Geometry for display.C
106
107     const Int_t kColorArm1   = kBlue ;
108
109     AliEMCALGeometry * geom = GetGeometry();
110
111     TString gn(geom->GetName());
112     gn.ToUpper(); 
113
114     // Define the shape of the Calorimeter 
115     TNode * top = gAlice->GetGeometry()->GetNode("alice") ; // See AliceGeom/Nodes
116     TNode * envelopNode = 0;
117     char *envn = "Envelop1";
118     if(!gn.Contains("SHISH")){
119       new TTUBS(envn, "Tubs that contains arm 1", "void", 
120               geom->GetEnvelop(0) -10, // rmin 
121               geom->GetEnvelop(1) +40 ,// rmax
122               geom->GetEnvelop(2)/2.0, // half length in Z
123               geom->GetArm1PhiMin(),   // minimum phi angle
124               geom->GetArm1PhiMax()    // maximum phi angle
125         );
126       top->cd();
127       envelopNode = new TNode(envn, "Arm1 Envelop", "Envelop1", 0., 0., 0., "") ;
128     } else {
129       if(gn.Contains("WSUC")) {
130         envelopNode = BuildGeometryOfWSUC();
131       } else { // Shish-kebab now for compact and TRD1 cases (ALIC)
132         envn="Envelop2";
133         TPGON *pgon = new TPGON(envn, "PGON that contains arm 1", "void", 
134         geom->GetArm1PhiMin(),geom->GetArm1PhiMax()-geom->GetArm1PhiMin(),geom->GetNPhiSuperModule(), 2);
135       // define section
136         if (fEnvelop1.GetSize()<8) {
137           fEnvelop1.Set(10);
138           fEnvelop1[0] = geom->GetArm1PhiMin(); // minimum phi angle
139           fEnvelop1[1] = geom->GetArm1PhiMax() - geom->GetArm1PhiMin(); // angular range in phi
140           fEnvelop1[2] = geom->GetNPhi(); // number of sections in phi
141           fEnvelop1[3] = 2; // 2 z coordinates
142           fEnvelop1[4] = -350.;
143           fEnvelop1[5] = geom->GetEnvelop(0) ; // rmin at z1
144           fEnvelop1[6] = geom->GetEnvelop(1) ; // rmax at z1
145           fEnvelop1[7] = 350.;
146           fEnvelop1[8] = fEnvelop1[5] ;        // radii are the same.
147           fEnvelop1[9] = fEnvelop1[6] ;        // radii are the same.
148           
149           if(gn.Contains("SHISH"))
150             fEnvelop1[2] = geom->GetNPhiSuperModule();
151         }
152         pgon->DefineSection(0, fEnvelop1[4],  fEnvelop1[5], fEnvelop1[6]);
153         pgon->DefineSection(1, fEnvelop1[7],  fEnvelop1[5], fEnvelop1[6]);
154         top->cd();
155         envelopNode = new TNode(envn, "Arm1 Envelop2", envn, 0., 0., 0., "") ;
156       }
157     }
158                                 
159     envelopNode->SetLineColor(kColorArm1) ;
160     fNodes->Add(envelopNode);
161 }
162
163 //______________________________________________________________________
164 TNode *AliEMCALv0::BuildGeometryOfWSUC()
165
166   // June 8, 2005; see directory geant3/TGeant3/G3toRoot.cxx
167   // enum EColor { kWhite, kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan } - see $ROOTSYS/include/Gtypes.h
168    AliEMCALGeometry * g = GetGeometry(); 
169    TNode * top = gAlice->GetGeometry()->GetNode("alice") ; // See AliceGeom/Nodes
170    top->cd();
171
172    TNode *envelopNode = 0;
173    char *name = "";
174    /*
175     name = "WSUC";
176    new TBRIK(name, "WSUC(XEN1 in Geant)","void",fEnvelop1[0],fEnvelop1[1],fEnvelop1[2]);
177     envelopNode = new TNode(name, "envelope for WSUC", name, 0., 0., 0., "");
178    envelopNode->SetVisibility(0);
179    */
180
181    TNode *emod=0, *scmx=0;
182    name = "SMOD"; // super module
183    new TBRIK(name, "SMOD(SMOD in Geant)","void", fSmodPar0,fSmodPar1,fSmodPar2);
184    if(envelopNode) envelopNode->cd();
185    TNode *smod = new TNode(name, "SMOD", name, 0., 0., 0., "");   
186    smod->SetLineColor(kBlue) ;
187    if(envelopNode==0) envelopNode = smod;
188
189    name = "EMOD"; // see CreateEMOD 
190    TTRD1 *emodTrd1 = new TTRD1(name, "EMOD(EMOD in Geant)","void", float(fParEMOD[0]),
191    float(fParEMOD[1]),float(fParEMOD[2]),float(fParEMOD[3]));
192
193    // SCMX = EMOD/4 for simplicity of drawing
194    name = "SCMX";
195    Float_t dz=0.,theta=0.,phi=0.,h1=0.,bl1=0.,tl1=0.,alpha1=0.,h2=0.,bl2=0.,tl2=0.,alpha2=0.;
196    h1     = emodTrd1->GetDy()/2.;
197    bl1    = emodTrd1->GetDx()/2.;
198    tl1    = bl1;
199    alpha1 = 0.;
200    h2     = emodTrd1->GetDy()/2.;
201    bl2    = emodTrd1->GetDx2()/2.;
202    tl2    = bl2;
203    alpha2 = 0.;
204
205    dz       = emodTrd1->GetDz();
206    double dr = TMath::Sqrt((h2-h1)*(h2-h1)+(bl2-bl1)*(bl2-bl1));
207    theta    = TMath::ATan2(dr,2.*dz) * TMath::RadToDeg(); 
208    phi      = 180.;
209
210    TTRAP *scmxTrap = new TTRAP(name, "SCMX(SCMX as in Geant)","void",
211              dz,theta,phi, h1,bl1,tl1,alpha1, h2,bl2,tl2,alpha2);
212    //   scmxTrap->Dump();
213    Float_t xShiftSCMX = (emodTrd1->GetDx() + emodTrd1->GetDx2())/4.;
214    Float_t yShiftSCMX = emodTrd1->GetDy()/2.;
215    printf(" xShiftSCMX %7.4f yShiftSCMX %7.4f \n",xShiftSCMX,yShiftSCMX);
216
217    name = "EMOD"; // see CreateEMOD 
218    smod->cd();
219    
220    AliEMCALShishKebabTrd1Module *mod=0;
221    Double_t  angle=90., xpos=0.,ypos=0.,zpos=0., xposSCMX=0.,yposSCMX=0.,zposSCMX=0.;
222    char rtmn[100], rtmt[100];
223    TRotMatrix *rtm=0, *rtmSCMX=0;
224    int numEmod=0;
225    for(int iz=0; iz<g->GetNZ(); iz++) {
226      mod   = (AliEMCALShishKebabTrd1Module*)fShishKebabModules->At(iz);
227      zpos = mod->GetPosZ()      - fSmodPar2;
228      ypos = mod->GetPosXfromR() - fSmodPar1;
229
230      angle = mod->GetThetaInDegree();
231      sprintf(rtmn,"rmEmod%5.1f",angle);
232      sprintf(rtmt,"rotation matrix for EMOD, iz=%i, angle = %6.3f",iz, angle);
233      if(iz==0) rtm = new TRotMatrix(rtmn, rtmt,0.,0., 90.,0., 90.,90.); // z'(x); y'(y); x'(z)
234      else      rtm = new TRotMatrix(rtmn, rtmt,90.-angle,270., 90.0,0.0, angle,90.);
235
236      TGeometry *tg = gAlice->GetGeometry();
237      for(int ix=0; ix<g->GetNPhi(); ix++) { // flat in phi
238        xpos = g->GetPhiModuleSize()*(2*ix+1 - g->GetNPhi())/2.;
239        sprintf(rtmt,"EMOD, iz %i, ix %i, angle %6.3f",iz,ix, angle);
240        TString namNode=name;
241        namNode += numEmod++;
242        smod->cd();
243        emod = new TNode(namNode.Data(), rtmt, (TShape*)emodTrd1, xpos,ypos,zpos,rtm);   
244        //       emod->SetLineColor(kGreen) ;
245        emod->SetVisibility(0); // SCMX will bi visible 
246        if(scmxTrap) { // 4(2x2) sensetive volume inside EMOD
247          emod->cd();
248          zposSCMX = 0.;
249          for(int jy=0; jy<2; jy++){ // division on y
250            yposSCMX = yShiftSCMX *(2*jy - 1); 
251            for(int jx=0; jx<2; jx++){ // division on x
252              Double_t theta1=90.,phi1=0., theta2=90.,phi2=90., theta3=0.,phi3=0 ; 
253              xposSCMX = xShiftSCMX *(2*jx - 1);              
254              namNode = "SCMX";
255              namNode += jy;
256              namNode += jx;
257              sprintf(rtmn,"rm%s",namNode.Data());
258              sprintf(rtmt,"rotation matrix for %s inside EMOD",namNode.Data());
259              rtmSCMX = tg->GetRotMatrix(rtmn);
260              if(jx == 1) {
261                phi1 = 180.; // x' = -x
262                phi2 = 270.; // y' = -y
263              }
264              if(rtmSCMX == 0) rtmSCMX = new TRotMatrix(rtmn,rtmt, theta1,phi1, theta2,phi2, theta3,phi3);
265              sprintf(rtmt,"%s inside %s", namNode.Data(), emod->GetName());
266              scmx = new TNode(namNode.Data(), rtmt, (TShape*)scmxTrap, xposSCMX,yposSCMX,zposSCMX,rtmSCMX);   
267              scmx->SetLineColor(kMagenta);
268            }
269          }
270        }
271      }
272    }
273    //   emod->Draw(); // for testing
274
275    return envelopNode;
276 }
277
278 //______________________________________________________________________
279 void AliEMCALv0::CreateGeometry()
280 {
281   // Create the EMCAL geometry for Geant
282   // Geometry of a tower
283   //|-----------------------------------------------------| XEN1
284   //| |                                                 | |
285   //| |    Al thickness = GetAlFrontThickness()         | |
286   //| |                                                 | |
287   //| |                                                 | |
288   //| |                                                 | |
289   //|  -------------------------------------------------  |
290   //| |    Air Gap = GetGap2Active()                    | |
291   //| |                                                 | |
292   //|  -------------------------------------------------  |
293   //| |    XU1 : XPST (ECAL e = GetECScintThick() )     | |
294   //|  -------------------------------------------------  |
295   //| |    XU1 : XPBX (ECAL e = GetECPbRadThick() )     | |
296   //|  -------------------------------------------------  |
297   //| |    XU1 : XPST (ECAL e = GetECScintThick()       | |
298   //|  -------------------------------------------------  |
299   //| |    XU1 : XPBX (ECAL e = GetECPbRadThick() )     | |
300   //|  -------------------------------------------------  |
301   //|    etc ..... GetNECLayers() - 1 times               |
302   //|  -------------------------------------------------  |
303   //| |    XUNLayer : XPST (ECAL e = GetECScintThick()  | |
304   //|  -------------------------------------------------  |
305
306     AliEMCALGeometry * geom = GetGeometry() ; 
307     TString gn(geom->GetName());
308     gn.ToUpper(); 
309
310     if(!(geom->IsInitialized())){
311         Error("CreateGeometry","EMCAL Geometry class has not been set up.");
312     } // end if
313
314     // Get pointer to the array containing media indices
315     fIdTmedArr = fIdtmed->GetArray() - 1599 ;
316
317     fIdRotm = 1;
318     //  gMC->Matrix(nmat, theta1, phi1, theta2, phi2, theta3, phi3) - see AliModule
319     AliMatrix(fIdRotm, 90.0, 0., 90.0, 90.0, 0.0, 0.0) ; 
320
321     // Create the EMCAL Mother Volume (a polygone) within which to place the Detector and named XEN1 
322
323     Float_t envelopA[10];
324     if(gn.Contains("TRD1") && gn.Contains("WSUC") ) { // TRD1 for WSUC facility
325       // 17-may-05 - just BOX
326       envelopA[0] = 26;
327       envelopA[1] = 15;
328       envelopA[2] = 30;
329       gMC->Gsvolu("XEN1", "BOX", fIdTmedArr[kIdSC], envelopA, 3) ;
330       fEnvelop1.Set(3);
331       for(int i=0; i<3; i++) fEnvelop1[i] = envelopA[i]; // 23-may-05  
332       // Position the EMCAL Mother Volume (XEN1) in WSUC  
333       gMC->Gspos("XEN1", 1, "WSUC", 0.0, 0.0, 0.0, fIdRotm, "ONLY") ;
334     } else { 
335       envelopA[0] = geom->GetArm1PhiMin();                         // minimum phi angle
336       envelopA[1] = geom->GetArm1PhiMax() - geom->GetArm1PhiMin(); // angular range in phi
337       envelopA[2] = geom->GetNPhi();                               // number of sections in phi
338       envelopA[3] = 2;                                             // 2 z coordinates
339       // envelopA[4] = geom->ZFromEtaR(geom->GetEnvelop(1),
340 //geom->GetArm1EtaMin());       // z coordinate 1
341       envelopA[4] = -350.; // AM 30/5/2006
342       
343     //add some padding for mother volume
344       envelopA[5] = geom->GetEnvelop(0) ;                          // rmin at z1
345       envelopA[6] = geom->GetEnvelop(1) ;                          // rmax at z1
346 //      envelopA[7] = geom->ZFromEtaR(geom->GetEnvelop(1),
347 //                                geom->GetArm1EtaMax());        // z coordinate 2
348       envelopA[7] = 350.; // AM 30/5/2006
349       
350       envelopA[8] = envelopA[5] ;                                  // radii are the same.
351       envelopA[9] = envelopA[6] ;                                  // radii are the same.
352
353       if(gn.Contains("SHISH")) envelopA[2] = geom->GetNPhiSuperModule();
354
355       gMC->Gsvolu("XEN1", "PGON", fIdTmedArr[kIdAIR], envelopA, 10) ;   // Polygone filled with air 
356       fEnvelop1.Set(10, envelopA);
357       if (gDebug==2) {
358         printf("CreateGeometry: XEN1 = %f, %f\n", envelopA[5], envelopA[6]); 
359         printf("CreateGeometry: XU0 = %f, %f\n", envelopA[5], envelopA[6]); 
360       }
361     // Position the EMCAL Mother Volume (XEN1) in Alice (ALIC)  
362       gMC->Gspos(geom->GetNameOfEMCALEnvelope(), 1, "ALIC", 0.0, 0.0, 0.0, fIdRotm, "ONLY") ;
363     }
364
365     if(gn.Contains("SHISH")){
366       // COMPACT, TRD1
367       AliDebug(2,Form("Shish-Kebab geometry : %s", GetTitle())); 
368       CreateShishKebabGeometry();
369     }
370 }
371
372 //______________________________________________________________________
373 void AliEMCALv0::Init(void)
374 {
375     // Just prints an information message
376   AliEMCAL::Init();
377   if(AliLog::GetGlobalDebugLevel()>0) { 
378     TString message("\n") ; 
379     message += "*****************************************\n" ;
380     
381     // Here the EMCAL initialisation code (if any!)
382     
383     AliEMCALGeometry * geom = GetGeometry() ; 
384     
385     if (geom!=0) {   
386       message += "AliEMCAL " ; 
387       message += Version() ; 
388       message += "EMCAL geometry initialized for " ; 
389       message += geom->GetName()  ;
390     }
391     else {
392       message += "AliEMCAL " ; 
393       message += Version() ;  
394       message += "EMCAL geometry initialization failed !" ; 
395     }
396     message += "\n*****************************************" ;
397     printf(message.Data() ) ; 
398   }
399 }
400
401 // 24-aug-04 by PAI
402 //______________________________________________________________________
403 void AliEMCALv0::CreateShishKebabGeometry()
404 {  
405   // TRD1
406   AliEMCALGeometry * g = GetGeometry(); 
407   TString gn(g->GetName()); gn.ToUpper(); 
408   // see AliModule::fFIdTmedArr
409   //  fIdTmedArr = fIdtmed->GetArray() - 1599 ; // see AliEMCAL::::CreateMaterials()
410   //  int kIdAIR=1599, kIdPB = 1600, kIdSC = 1601, kIdSTEEL = 1603;
411   //  idAL = 1602;
412   Double_t par[10], xpos=0., ypos=0., zpos=0.;
413
414   CreateSmod(g->GetNameOfEMCALEnvelope());
415
416   CreateEmod("SMOD","EMOD"); // 18-may-05
417
418   if(g->GetKey110DEG()) CreateEmod("SM10","EMOD"); // Nov 1,2006 
419
420   // Sensitive SC  (2x2 tiles)
421   double parSCM0[5]={0,0,0,0}, *dummy = 0, parTRAP[11];
422   Double_t trd1Angle = g->GetTrd1Angle()*TMath::DegToRad();
423    if(!gn.Contains("TRD")) { // standard module
424     par[0] = (g->GetECPbRadThick()+g->GetECScintThick())*g->GetNECLayers()/2.;
425     par[1] = par[2] = g->GetPhiTileSize();   // Symetric case
426     gMC->Gsvolu("SCM0", "BOX", fIdTmedArr[kIdSC], par, 3); // 2x2 tiles
427     gMC->Gspos("SCM0", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
428     // Division to tile size
429     gMC->Gsdvn("SCM1","SCM0", g->GetNPHIdiv(), 2); // y-axis
430     gMC->Gsdvn("SCM2","SCM1", g->GetNETAdiv(), 3); // z-axis
431   // put LED to the SCM2 
432     par[0] = g->GetECPbRadThick()/2.;
433     par[1] = par[2] = g->GetPhiTileSize()/2.; // Symetric case
434     gMC->Gsvolu("PBTI", "BOX", fIdTmedArr[kIdPB], par, 3);
435
436     printf(" Pb tiles \n");
437     int nr=0;
438     ypos = zpos = 0.0;
439     xpos = -fSampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
440     for(int ix=0; ix<g->GetNECLayers(); ix++){
441       gMC->Gspos("PBTI", ++nr, "SCM2", xpos, ypos, zpos, 0, "ONLY") ;
442     //    printf(" %i xpos %f \n", ix+1, xpos);
443       xpos += fSampleWidth;
444     } 
445     printf(" Number of Pb tiles in SCM2 %i \n", nr);
446   } else if(gn.Contains("TRD1")) { // TRD1 - 30-sep-04
447     //      double wallThickness = g->GetPhiModuleSize()/2. -  g->GetPhiTileSize(); // Need check
448     double wallThickness = g->GetPhiModuleSize()/g->GetNPHIdiv() -  g->GetPhiTileSize();
449     for(int i=0; i<3; i++) parSCM0[i] = fParEMOD[i] - wallThickness;
450     parSCM0[3] = fParEMOD[3];
451     gMC->Gsvolu("SCM0", "TRD1", fIdTmedArr[kIdAIR], parSCM0, 4);
452     gMC->Gspos("SCM0", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
453
454     if(g->GetNPHIdiv()==2 && g->GetNETAdiv()==2) {
455     // Division to tile size - 1-oct-04
456       AliDebug(2,Form(" Divide SCM0 on y-axis %i\n", g->GetNETAdiv()));
457       gMC->Gsdvn("SCMY","SCM0", g->GetNETAdiv(), 2); // y-axis
458     // Trapesoid 2x2
459       parTRAP[0] = parSCM0[3];    // dz
460       parTRAP[1] = TMath::ATan2((parSCM0[1]-parSCM0[0])/2.,2.*parSCM0[3])*180./TMath::Pi(); // theta
461       parTRAP[2] = 0.;           // phi
462       // bottom
463       parTRAP[3] = parSCM0[2]/2.; // H1
464       parTRAP[4] = parSCM0[0]/2.; // BL1
465       parTRAP[5] = parTRAP[4];    // TL1
466       parTRAP[6] = 0.0;           // ALP1
467       // top
468       parTRAP[7] = parSCM0[2]/2.; // H2
469       parTRAP[8] = parSCM0[1]/2.; // BL2
470       parTRAP[9] = parTRAP[8];    // TL2
471       parTRAP[10]= 0.0;           // ALP2
472       AliDebug(2,Form(" ** TRAP ** \n"));
473       for(int i=0; i<11; i++) AliDebug(3, Form(" par[%2.2i] %9.4f\n", i, parTRAP[i]));
474
475       gMC->Gsvolu("SCMX", "TRAP", fIdTmedArr[kIdSC], parTRAP, 11);
476       xpos = +(parSCM0[1]+parSCM0[0])/4.;
477       gMC->Gspos("SCMX", 1, "SCMY", xpos, 0.0, 0.0, 0, "ONLY") ;
478
479       // Using rotation because SCMX should be the same due to Pb tiles
480       xpos = -xpos; 
481       AliMatrix(fIdRotm, 90.0,180., 90.0, 270.0, 0.0,0.0) ;
482       gMC->Gspos("SCMX", 2, "SCMY", xpos, 0.0, 0.0, fIdRotm, "ONLY");
483     // put LED to the SCM0 
484       AliEMCALShishKebabTrd1Module *mod = (AliEMCALShishKebabTrd1Module*)fShishKebabModules->At(0);
485       gMC->Gsvolu("PBTI", "BOX", fIdTmedArr[kIdPB], dummy, 0);
486
487       par[1] = parSCM0[2]/2;            // y 
488       par[2] = g->GetECPbRadThick()/2.; // z
489
490       int nr=0;
491       ypos = 0.0; 
492       zpos = -fSampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
493       double xCenterSCMX =  (parTRAP[4] +  parTRAP[8])/2.;
494       if(!gn.Contains("NOPB")) { // for testing - 11-jul-05
495         AliDebug(2,Form(" Pb tiles \n"));
496         for(int iz=0; iz<g->GetNECLayers(); iz++){
497           par[0] = (parSCM0[0] + mod->GetTanBetta()*fSampleWidth*iz)/2.;
498           xpos   = par[0] - xCenterSCMX;
499           gMC->Gsposp("PBTI", ++nr, "SCMX", xpos, ypos, zpos, 0, "ONLY", par, 3) ;
500           AliDebug(3,Form(" %i xpos %f zpos %f par[0] %f \n", iz+1, xpos, zpos, par[0]));
501           zpos += fSampleWidth;
502         } 
503         AliDebug(2,Form(" Number of Pb tiles in SCMX %i \n", nr));
504       }
505     } else if(g->GetNPHIdiv()==3 && g->GetNETAdiv()==3) {
506       printf(" before AliEMCALv0::Trd1Tower3X3() : parSCM0");
507       for(int i=0; i<4; i++) printf(" %7.4f ", parSCM0[i]);
508       printf("\n"); 
509       Trd1Tower3X3(parSCM0);
510     } else if(g->GetNPHIdiv()==1 && g->GetNETAdiv()==1) {
511       // no division in SCM0
512       Trd1Tower1X1(parSCM0);
513     } else if(g->GetNPHIdiv()==4 && g->GetNETAdiv()==4) {
514       Trd1Tower4X4();
515     }
516   }
517 }
518
519 //______________________________________________________________________
520 void AliEMCALv0::CreateSmod(const char* mother)
521
522   // 18-may-05; mother="XEN1"; 
523   // child="SMOD" from first to 10th, "SM10" (11th and 12th) (TRD1 case)
524   AliEMCALGeometry * g = GetGeometry(); 
525   TString gn(g->GetName()); gn.ToUpper();
526
527   Double_t par[3], xpos=0., ypos=0., zpos=0., rpos=0., dphi=0., phi=0.0, phiRad=0.;
528   Double_t par1C = 0.;
529   //  ===== define Super Module from air - 14x30 module ==== ;
530   fSampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
531   AliDebug(2,Form("\n ## Super Module | fSampleWidth %5.3f ## %s \n", fSampleWidth, gn.Data()));
532   par[0] = g->GetShellThickness()/2.;
533   par[1] = g->GetPhiModuleSize()*g->GetNPhi()/2.; 
534   par[2] = g->GetEtaModuleSize()*15.; 
535   fIdRotm=0;
536   int nphism = g->GetNumberOfSuperModules()/2; // 20-may-05
537   if(nphism>0) {
538     dphi = (g->GetArm1PhiMax() - g->GetArm1PhiMin())/nphism;
539     //    if(g->GetKey110DEG()) dphi = (g->GetArm1PhiMax() - g->GetArm1PhiMin())/(nphism-1);
540     rpos = (g->GetEnvelop(0) + g->GetEnvelop(1))/2.;
541     AliDebug(2,Form(" rpos %8.2f : dphi %6.1f degree \n", rpos, dphi));
542   }
543
544   if(gn.Contains("WSUC")) {
545     par[0] = g->GetPhiModuleSize()*g->GetNPhi()/2.; 
546     par[1] = g->GetShellThickness()/2.;
547     par[2] = g->GetEtaModuleSize()*g->GetNZ()/2. + 5; 
548
549     gMC->Gsvolu("SMOD", "BOX", fIdTmedArr[kIdAIR], par, 3);
550
551     AliDebug(2,Form("SMOD in WSUC : tmed %i | dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", 
552                     fIdTmedArr[kIdAIR], par[0],par[1],par[2]));
553     fSmodPar0 = par[0]; 
554     fSmodPar1 = par[1];
555     fSmodPar2 = par[2];
556     nphism   =  g->GetNumberOfSuperModules();
557   } else {
558     if(gn.Contains("TRD")) {
559       par[2]  = 350./2.; // 11-oct-04 - for 26 division
560       if(gn.Contains("TRD1")) {
561         AliDebug(2,Form(" par[0] %7.2f (old) \n",  par[0]));
562         Float_t *parSM = g->GetSuperModulesPars(); 
563         for(int i=0; i<3; i++) par[i] = parSM[i];
564       }
565     }
566     gMC->Gsvolu("SMOD", "BOX", fIdTmedArr[kIdAIR], par, 3);
567     AliDebug(2,Form("tmed %i | dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", 
568                     fIdTmedArr[kIdAIR], par[0],par[1],par[2]));
569     fSmodPar0 = par[0]; 
570     fSmodPar2 = par[2];
571     if(g->GetKey110DEG()) { // 12-oct-05
572       par1C = par[1];
573       par[1] /= 2.;
574       gMC->Gsvolu("SM10", "BOX", fIdTmedArr[kIdAIR], par, 3);
575       AliDebug(2,Form(" Super module with name \"SM10\" was created too par[1] = %f\n", par[1]));
576       par[1] = par1C;
577     }
578   // Steel plate
579     if(g->GetSteelFrontThickness() > 0.0) { // 28-mar-05
580       par[0] = g->GetSteelFrontThickness()/2.;
581       gMC->Gsvolu("STPL", "BOX", fIdTmedArr[kIdSTEEL], par, 3);
582       printf("tmed %i | dx %7.2f dy %7.2f dz %7.2f (STPL) \n", fIdTmedArr[kIdSTEEL], par[0],par[1],par[2]);
583       xpos = -(g->GetShellThickness() - g->GetSteelFrontThickness())/2.;
584       gMC->Gspos("STPL", 1, "SMOD", xpos, 0.0, 0.0, 0, "ONLY") ;
585     }
586   }
587
588   int nr=0, nrsmod=0, i0=0;
589
590   // Turn whole super module
591   for(int i=i0; i<nphism; i++) {
592     if(gn.Contains("WSUC")) {
593       xpos = ypos = zpos = 0.0;
594       fIdRotm = 0;
595       gMC->Gspos("SMOD", 1, mother, xpos, ypos, zpos, fIdRotm, "ONLY") ;
596       printf(" fIdRotm %3i phi %6.1f(%5.3f) xpos %7.2f ypos %7.2f zpos %7.2f \n", 
597       fIdRotm, phi, phiRad, xpos, ypos, zpos);
598       nr++;
599     } else { // TRD1 
600       TString smName("SMOD"); // 12-oct-05
601       if(i==5 && g->GetKey110DEG()) {
602         smName = "SM10";
603         nrsmod = nr;
604         nr     = 0;
605       }
606       phi    = g->GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 70, 90, 110, 130, 150, 170
607       phiRad = phi*TMath::Pi()/180.;
608
609       AliMatrix(fIdRotm, 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
610
611       xpos = rpos * TMath::Cos(phiRad);
612       ypos = rpos * TMath::Sin(phiRad);
613       zpos = fSmodPar2; // 21-sep-04
614       if(i==5 && g->GetKey110DEG()) {
615         xpos += (par1C/2. * TMath::Sin(phiRad)); 
616         ypos -= (par1C/2. * TMath::Cos(phiRad)); 
617       }
618       
619       // 1th module in z-direction;
620       gMC->Gspos(smName.Data(), ++nr, mother, xpos, ypos, zpos, fIdRotm, "ONLY") ;
621       AliDebug(3, Form(" %s : %2i fIdRotm %3i phi %6.1f(%5.3f) xpos %7.2f ypos %7.2f zpos %7.2f : i %i \n", 
622                        smName.Data(), nr, fIdRotm, phi, phiRad, xpos, ypos, zpos, i));
623       // 2th module in z-direction;
624       if(gn.Contains("TRD")) {
625       // turn arround X axis; 0<phi<360
626         double phiy = 90. + phi + 180.;
627         if(phiy>=360.) phiy -= 360.;
628  
629         AliMatrix(fIdRotm, 90.0, phi, 90.0, phiy, 180.0, 0.0);
630         gMC->Gspos(smName.Data(), ++nr, mother, xpos, ypos, -zpos, fIdRotm, "ONLY");
631         AliDebug(3, Form(" %s : %2i fIdRotm %3i phiy %6.1f  xpos %7.2f ypos %7.2f zpos %7.2f \n", 
632                          smName.Data(), nr, fIdRotm, phiy, xpos, ypos, -zpos));
633       } else {
634         gMC->Gspos("SMOD", ++nr, mother, xpos, ypos, -zpos, fIdRotm, "ONLY");
635       }
636     }
637   }
638   AliDebug(2,Form(" Number of Super Modules %i \n", nr+nrsmod));
639 }
640
641 //______________________________________________________________________
642 void AliEMCALv0::CreateEmod(const char* mother, const char* child)
643
644   // 17-may-05; mother="SMOD"; child="EMOD"
645   AliEMCALGeometry * g = GetGeometry(); 
646   TString gn(g->GetName()); gn.ToUpper(); 
647   // Module definition
648   Double_t par[10], xpos=0., ypos=0., zpos=0.;
649   Double_t parSCPA[5], zposSCPA=0.; // passive SC - 13-MAY-05, TRD1 case
650   Double_t trd1Angle = g->GetTrd1Angle()*TMath::DegToRad(), tanTrd1 = TMath::Tan(trd1Angle/2.);
651   int nr=0;
652   fIdRotm=0;
653   if(!gn.Contains("TRD")) { // standard module
654     par[0] = (fSampleWidth*g->GetNECLayers())/2.; 
655     par[1] = par[2] = g->GetPhiModuleSize()/2.;
656     gMC->Gsvolu(child, "BOX", fIdTmedArr[kIdSTEEL], par, 3);
657
658   } else if (gn.Contains("TRD1")){ // TRD1 system coordinate iz differnet
659     if(strcmp(mother,"SMOD")==0) {
660       fParEMOD[0] = g->GetEtaModuleSize()/2.;   // dx1
661       fParEMOD[1] = g->Get2Trd1Dx2()/2.;        // dx2
662       fParEMOD[2] = g->GetPhiModuleSize()/2.;;  // dy
663       fParEMOD[3] = g->GetLongModuleSize()/2.;  // dz
664       gMC->Gsvolu(child, "TRD1", fIdTmedArr[kIdSTEEL], fParEMOD, 4);
665       //      if(gn.Contains("WSUC") || gn.Contains("MAY05")){
666       if(0){ // Jul 12 - should be checked
667         parSCPA[0] = g->GetEtaModuleSize()/2. + tanTrd1*g->GetFrontSteelStrip();   // dx1
668         parSCPA[1] = parSCPA[0]               + tanTrd1*g->GetPassiveScintThick(); // dx2
669         parSCPA[2] = g->GetPhiModuleSize()/2.;     // dy
670         parSCPA[3] = g->GetPassiveScintThick()/2.; // dz
671         gMC->Gsvolu("SCPA", "TRD1", fIdTmedArr[kIdSC], parSCPA, 4);
672         zposSCPA   = -fParEMOD[3] + g->GetFrontSteelStrip() + g->GetPassiveScintThick()/2.;
673         gMC->Gspos ("SCPA", ++nr, child, 0.0, 0.0, zposSCPA, 0, "ONLY");
674       }
675     }
676   }
677
678   nr   = 0;
679   if(gn.Contains("TRD")) { // 30-sep-04; 27-jan-05 - as for TRD1
680     // X->Z(0, 0); Y->Y(90, 90); Z->X(90, 0)
681     AliEMCALShishKebabTrd1Module *mod=0; // current module
682
683     for(int iz=0; iz<g->GetNZ(); iz++) {
684       Double_t  angle=90., phiOK=0;
685       if(gn.Contains("TRD1")) {
686         mod = (AliEMCALShishKebabTrd1Module*)fShishKebabModules->At(iz);
687         angle = mod->GetThetaInDegree();
688         if(!gn.Contains("WSUC")) { // ALICE 
689           if(iz==0) AliMatrix(fIdRotm, 0.,0., 90.,90., 90.,0.); // z'(x); y'(y); x'(z)
690           else      AliMatrix(fIdRotm, 90.-angle,180., 90.0,90.0, angle, 0.);
691           phiOK = mod->GetCenterOfModule().Phi()*180./TMath::Pi(); 
692           //          printf(" %2i | angle | %6.3f - %6.3f = %6.3f(eta %5.3f)\n", 
693           //iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule());
694           xpos = mod->GetPosXfromR() + g->GetSteelFrontThickness() - fSmodPar0;
695           zpos = mod->GetPosZ() - fSmodPar2;
696
697           int iyMax = g->GetNPhi();
698           if(strcmp(mother,"SMOD") && g->GetKey110DEG()) {
699             iyMax /= 2;
700           }
701           for(int iy=0; iy<iyMax; iy++) { // flat in phi
702             ypos = g->GetPhiModuleSize()*(2*iy+1 - iyMax)/2.;
703             gMC->Gspos(child, ++nr, mother, xpos, ypos, zpos, fIdRotm, "ONLY") ;
704         //printf(" %2i xpos %7.2f ypos %7.2f zpos %7.2f fIdRotm %i\n", nr, xpos, ypos, zpos, fIdRotm);
705             AliDebug(3,Form("%3.3i(%2.2i,%2.2i) ", nr,iy+1,iz+1));
706           }
707           //PH          printf("\n");
708         } else {
709           if(iz==0) AliMatrix(fIdRotm, 0.,0., 90.,0., 90.,90.); // (x')z; y'(x); z'(y)
710           else      AliMatrix(fIdRotm, 90-angle,270., 90.0,0.0, angle,90.);
711           phiOK = mod->GetCenterOfModule().Phi()*180./TMath::Pi(); 
712           printf(" %2i | angle -phiOK | %6.3f - %6.3f = %6.3f(eta %5.3f)\n", 
713           iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule());
714           zpos = mod->GetPosZ()      - fSmodPar2;
715           ypos = mod->GetPosXfromR() - fSmodPar1;
716           printf(" zpos %7.2f ypos %7.2f fIdRotm %i\n xpos ", zpos, xpos, fIdRotm);
717           for(int ix=0; ix<g->GetNPhi(); ix++) { // flat in phi
718             xpos = g->GetPhiModuleSize()*(2*ix+1 - g->GetNPhi())/2.;
719             gMC->Gspos(child, ++nr, mother, xpos, ypos, zpos, fIdRotm, "ONLY") ;
720             printf(" %7.2f ", xpos);
721           }
722           printf("\n");
723         }
724       }
725     } 
726   } else {
727     xpos = g->GetSteelFrontThickness()/2.;
728     for(int iz=0; iz<g->GetNZ(); iz++) {
729       zpos = -fSmodPar2 + g->GetEtaModuleSize()*(2*iz+1)/2.;
730       for(int iy=0; iy<g->GetNPhi(); iy++) {
731         ypos = g->GetPhiModuleSize()*(2*iy+1 - g->GetNPhi())/2.;
732         gMC->Gspos(child, ++nr, mother, xpos, ypos, zpos, 0, "ONLY") ;
733       //printf(" %2i xpos %7.2f ypos %7.2f zpos %7.2f \n", nr, xpos, ypos, zpos);
734       }
735     }
736   }
737   AliDebug(2,Form(" Number of modules in Super Module(%s) %i \n", mother, nr));
738 }
739
740 //______________________________________________________________________
741 void AliEMCALv0::Trd1Tower3X3(const double *parSCM0)
742 {
743   // Started Dec 8,2004 by PAI
744   // Fixed Nov 13,2006
745   printf(" AliEMCALv0::Trd1Tower3X3() : parSCM0");
746   for(int i=0; i<4; i++) printf(" %7.4f ", parSCM0[i]);
747   printf("\n"); 
748   // Nov 10, 2006 - different name of SCMX
749   double parTRAP[11], *dummy=0;
750   AliEMCALGeometry * g = GetGeometry(); 
751   TString gn(g->GetName()), scmx; 
752   gn.ToUpper(); 
753  // Division to tile size 
754   AliDebug(2,Form("Trd1Tower3X3() : Divide SCM0 on y-axis %i", g->GetNETAdiv()));
755   gMC->Gsdvn("SCMY","SCM0", g->GetNETAdiv(), 2); // y-axis
756   double dx1=parSCM0[0], dx2=parSCM0[1], dy=parSCM0[2], dz=parSCM0[3];
757   double ndiv=3., xpos=0.0;
758   // should be defined once
759   gMC->Gsvolu("PBTI", "BOX", fIdTmedArr[kIdPB], dummy, 0);
760   for(int ix=1; ix<=3; ix++) { // 3X3
761     scmx = "SCX"; // Nov 10,2006 
762     // ix=1
763     parTRAP[0] = dz;
764     double xCentBot = 2.*dx1/3.;
765     double xCentTop = 2.*(dx2/4. + dx1/12.);
766     parTRAP[1] = TMath::ATan2((xCentTop-xCentBot),2.*dz)*TMath::RadToDeg(); // theta
767     parTRAP[2] = 0.;           // phi
768     // bottom
769     parTRAP[3] = dy/ndiv;      // H1
770     parTRAP[4] = dx1/ndiv;     // BL1
771     parTRAP[5] = parTRAP[4];   // TL1
772     parTRAP[6] = 0.0;          // ALP1
773     // top
774     parTRAP[7] = dy/ndiv;      // H2
775     parTRAP[8] = dx2/2 - dx1/6.;// BL2
776     parTRAP[9] = parTRAP[8];   // TL2
777     parTRAP[10]= 0.0;          // ALP2
778     xpos = (xCentBot+xCentTop)/2.;
779
780     if      (ix==3) {
781       parTRAP[1] = -parTRAP[1];
782       xpos = -xpos;
783     } else if(ix==2) { // central part is box but we treat as trapesoid due to numbering
784       parTRAP[1] = 0.;
785       parTRAP[8] = dx1/ndiv;     // BL2
786       parTRAP[9] = parTRAP[8];   // TL2
787       xpos = 0.0;
788     }
789     AliDebug(2,Form(" ** TRAP ** xpos %9.3f\n", xpos));
790     for(int i=0; i<11; i++) AliDebug(2,Form(" par[%2.2i] %9.4f\n", i, parTRAP[i]));
791
792     scmx += ix;
793     gMC->Gsvolu(scmx.Data(), "TRAP", fIdTmedArr[kIdSC], parTRAP, 11);
794     gMC->Gspos(scmx.Data(), 1, "SCMY", xpos, 0.0, 0.0, 0, "ONLY") ;
795
796     PbInTrap(parTRAP, scmx);
797   }
798   AliDebug(2,Form("Trd1Tower3X3()", "Ver. 1.0 : was tested."));
799 }
800
801 // 8-dec-04 by PAI
802 //______________________________________________________________________
803 void AliEMCALv0::PbInTrap(const double parTRAP[11], TString n)
804 {
805  // see for example CreateShishKebabGeometry(); just for case TRD1
806   static int nr=0;
807   AliDebug(2,Form(" Pb tiles : nrstart %i\n", nr));
808   AliEMCALGeometry * g = GetGeometry(); 
809
810   double par[3];
811   //  double fSampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
812   double xpos = 0.0, ypos = 0.0;
813   double zpos = -fSampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
814
815   double coef = (parTRAP[8] -  parTRAP[4]) / (2.*parTRAP[0]);
816   double xCenterSCMX =  (parTRAP[4] +  parTRAP[8])/2.; // ??
817   //  double tan = TMath::Tan(parTRAP[1]*TMath::DegToRad());
818
819   par[1] = parTRAP[3];              // y 
820   par[2] = g->GetECPbRadThick()/2.; // z
821   for(int iz=0; iz<g->GetNECLayers(); iz++){
822     par[0] = parTRAP[4] + coef*fSampleWidth*iz;
823     xpos   = par[0] - xCenterSCMX;
824     if(parTRAP[1] < 0.) xpos = -xpos;
825     gMC->Gsposp("PBTI", ++nr, n.Data(), xpos, ypos, zpos, 0, "ONLY", par, 3) ;
826     AliDebug(2,Form(" %i xpos %9.3f zpos %9.3f par[0] %9.3f |", iz+1, xpos, zpos, par[0]));
827     zpos += fSampleWidth;
828     if(iz%2>0) printf("\n");
829   } 
830   AliDebug(2,Form(" Number of Pb tiles in SCMX %i coef %9.7f \n", nr, coef));
831   AliDebug(2,Form(" par[1] %9.3f  par[2] %9.3f ypos %9.3f \n", par[1], par[2], ypos)); 
832   AliDebug(2,Form(" PbInTrap Ver. 1.0 : was tested."));
833 }
834
835 // 8-dec-04 by PAI
836 //______________________________________________________________________
837 void AliEMCALv0::Trd1Tower4X4()
838 {
839  // Not ready yet
840 }
841
842 //______________________________________________________________________
843 void AliEMCALv0::Trd1Tower1X1(double *parSCM0)
844 {
845   // Started Nov 22,2006 by PAI
846   AliDebug(1," AliEMCALv0::Trd1Tower1X1() : parSCM0");
847   for(int i=0; i<4; i++) printf(" %7.4f ", parSCM0[i]);
848   printf("\n"); 
849
850   // No division - keeping the same volume logic 
851   // and as consequence the same abs is scheme
852   AliDebug(2,"Trd1Tower1X1() : Create SCMX(SCMY) as SCM0");
853
854   gMC->Gsvolu("SCMY", "TRD1", fIdTmedArr[kIdAIR], parSCM0, 4);
855   gMC->Gspos("SCMY", 1, "SCM0", 0.0, 0.0, 0.0, 0, "ONLY");
856   gMC->Gsvolu("SCMX", "TRD1", fIdTmedArr[kIdSC], parSCM0, 4);
857   gMC->Gspos("SCMX", 1, "SCMY", 0.0, 0.0, 0.0, 0, "ONLY");
858
859   // should be defined once
860   double *dummy=0;
861   gMC->Gsvolu("PBTI", "BOX", fIdTmedArr[kIdPB], dummy, 0);
862
863   PbInTrd1(parSCM0, "SCMX");
864
865   AliDebug(1,"Trd1Tower1X1() : Ver. 0.1 : was tested.");
866 }
867
868 //______________________________________________________________________
869 void AliEMCALv0::PbInTrd1(double *parTrd1, TString n)
870 {
871  // see PbInTrap(const double parTrd1[11], TString n)
872   static int nr=0, ndeb=2;
873   AliDebug(ndeb,Form(" Pb tiles : nrstart %i\n", nr));
874   AliEMCALGeometry * g = GetGeometry(); 
875
876   double par[3];
877   //  double fSampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
878   double xpos = 0.0, ypos = 0.0;
879   double zpos = -fSampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
880   double coef = (parTrd1[1] -  parTrd1[0]) / (2.*parTrd1[3]);
881
882   par[1] = parTrd1[2];              // y 
883   par[2] = g->GetECPbRadThick()/2.; // z
884
885   for(int iz=0; iz<g->GetNECLayers(); iz++){
886     par[0] = parTrd1[0] + coef*fSampleWidth*iz;
887     gMC->Gsposp("PBTI", ++nr, n.Data(), xpos, ypos, zpos, 0, "ONLY", par, 3) ;
888     AliDebug(2,Form(" %i xpos %9.3f zpos %9.3f par[0] %9.3f |", iz+1, xpos, zpos, par[0]));
889     zpos += fSampleWidth;
890     if(iz%2>0) printf("\n");
891   } 
892   AliDebug(ndeb,Form(" Number of Pb tiles in SCMX %i coef %9.7f ", nr, coef));
893   AliDebug(ndeb,Form(" PbInTrd1 Ver. 0.1 : was tested."));
894 }
895
896 // 3-feb-05
897 //______________________________________________________________________
898 void AliEMCALv0::Scm0InTrd2(const AliEMCALGeometry * g, const Double_t emodPar[5], Double_t parSCM0[5])
899 {
900   // Passive material inside the detector
901   double wallThickness = g->GetPhiModuleSize()/2. -  g->GetPhiTileSize(); //Need check
902   AliDebug(2,Form(" wall thickness %7.5f \n", wallThickness));
903   for(int i=0; i<4; i++) { // on pictures sometimes I can not see 0 -> be carefull!!
904     parSCM0[i] = emodPar[i] - wallThickness;
905     AliDebug(2,Form(" %i parSCMO %7.3f emodPar %7.3f : dif %7.3f \n", 
906                     i, parSCM0[i],emodPar[i], parSCM0[i]-emodPar[i]));
907   }
908   parSCM0[4] = emodPar[4];
909   gMC->Gsvolu("SCM0", "TRD2", fIdTmedArr[kIdSC], parSCM0, 5); // kIdAIR -> kIdSC
910   gMC->Gspos("SCM0", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
911   // Division 
912   if(g->GetNPHIdiv()==2 && g->GetNETAdiv()==2) {
913     Division2X2InScm0(g, parSCM0);
914   } else {
915     Info("Scm0InTrd2"," no division SCM0 in this geometry |%s|\n", g->GetName());
916     assert(0);
917   }
918 }
919
920 //______________________________________________________________________
921 void AliEMCALv0::Division2X2InScm0(const AliEMCALGeometry * g, const Double_t parSCM0[5])
922 {
923   // Division 2X2
924   Double_t parTRAP[11], xpos=0.,ypos=0., dx1=0.0,dx2=0.,dy1=0.0,dy2=0.,dz=0.
925   ,dr1=0.0,dr2=0.;
926   fIdRotm=0;
927
928   Info("Division2X2InScm0","Divide SCM0 on y-axis %i\n", g->GetNETAdiv());
929   TString n("SCMX"), overLapFlagSCMY("ONLY"), overLapFlagSCMX("ONLY");
930   n = "SCM0"; // for testing - 14-mar-05
931   if(n=="SCM0"){
932     PbInTrapForTrd2(parSCM0, n);
933     // overLapFlagSCMY=overLapFlagSCMX="MANY"; // do not work
934     return;
935   }
936
937   dy1 = parSCM0[2] , dy2 = parSCM0[3], dz = parSCM0[4];
938
939   parTRAP[0] = parSCM0[4];    // dz
940   parTRAP[1] = TMath::ATan2((dy2-dy1)/2.,2.*dz)*TMath::RadToDeg();
941   parTRAP[2] = 90.;           // phi
942   // bottom
943   parTRAP[3] = parSCM0[2]/2.; // H1
944   parTRAP[4] = parSCM0[0];    // BL1
945   parTRAP[5] = parTRAP[4];    // TL1
946   parTRAP[6] = 0.0;           // ALP1
947   // top
948   parTRAP[7] = parSCM0[3]/2.; // H2
949   parTRAP[8] = parSCM0[1];    // BL2
950   parTRAP[9] = parTRAP[8];    // TL2
951   parTRAP[10]= 0.0;           // ALP2
952   AliDebug(2,Form(" ** SCMY ** \n"));
953            for(int i=0; i<11; i++) AliDebug(2,Form(" par[%2.2i] %9.4f\n", i, parTRAP[i]));
954
955   fIdRotm=0;
956   gMC->Gsvolu("SCMY", "TRAP", fIdTmedArr[kIdSC], parTRAP, 11); // kIdAIR -> kIdSC
957   ypos = +(parTRAP[3]+parTRAP[7])/2.; //
958   AliDebug(2,Form(" Y shift SCMY inside SCM0 : %7.3f : opt %s\n", ypos, overLapFlagSCMY.Data())); 
959   gMC->Gspos("SCMY", 1, "SCM0", 0.0, ypos, 0.0, fIdRotm,  overLapFlagSCMY.Data()) ;
960   // Rotation SCMY around z-axis on 180 degree; x'=-x; y'=-y and z=z
961   AliMatrix(fIdRotm, 90.0,180., 90.0, 270.0, 0.0,0.0) ;
962   // We may have problem with numeration due to rotation - 4-feb-05
963   gMC->Gspos("SCMY", 2, "SCM0", 0.0, -ypos, 0.0, fIdRotm,  overLapFlagSCMY.Data()); 
964
965   Info("Division2X2InScm0","Divide SCMY on x-axis %i\n", g->GetNPHIdiv());
966   dx1 = parSCM0[0]; 
967   dx2 = parSCM0[1]; 
968   dr1=TMath::Sqrt(dx1*dx1+dy1*dy1);
969   dr2=TMath::Sqrt(dx2*dx2+dy2*dy2);
970
971   parTRAP[0] = parSCM0[4];    // dz
972   parTRAP[1] = TMath::ATan2((dr2-dr1)/2.,2.*dz)*TMath::RadToDeg(); // 
973   parTRAP[2] = 45.;           // phi
974   // bottom
975   parTRAP[3] = parSCM0[2]/2.; // H1
976   parTRAP[4] = parSCM0[0]/2.; // BL1
977   parTRAP[5] = parTRAP[4];    // TL1
978   parTRAP[6] = 0.0;           // ALP1
979   // top
980   parTRAP[7] = parSCM0[3]/2.; // H2
981   parTRAP[8] = parSCM0[1]/2;  // BL2
982   parTRAP[9] = parTRAP[8];    // TL2
983   parTRAP[10]= 0.0;           // ALP2
984   AliDebug(2,Form(" ** SCMX ** \n"));
985   for(int i=0; i<11; i++) AliDebug(2,Form(" par[%2.2i] %9.4f\n", i, parTRAP[i]));
986
987   fIdRotm=0;
988   gMC->Gsvolu("SCMX", "TRAP", fIdTmedArr[kIdSC], parTRAP, 11);
989   xpos = (parTRAP[4]+parTRAP[8])/2.;
990   AliDebug(2,Form(" X shift SCMX inside SCMX : %7.3f : opt %s\n", xpos, overLapFlagSCMX.Data())); 
991   gMC->Gspos("SCMX", 1, "SCMY", xpos, 0.0, 0.0, fIdRotm,  overLapFlagSCMX.Data()) ;
992   //  AliMatrix(fIdRotm, 90.0,270., 90.0, 0.0, 0.0,0.0); // x'=-y; y'=x; z'=z
993   AliMatrix(fIdRotm, 90.0,90., 90.0, -180.0, 0.0,0.0);     // x'=y;  y'=-x; z'=z
994   gMC->Gspos("SCMX", 2, "SCMY", -xpos, 0.0, 0.0, fIdRotm,  overLapFlagSCMX.Data()) ;
995   // PB:
996   if(n=="SCMX" && overLapFlagSCMY == "ONLY") {
997     PbInTrapForTrd2(parTRAP, n);
998   }
999
1000
1001 // 4-feb-05 by PAI
1002 //______________________________________________________________________
1003 void AliEMCALv0::PbInTrapForTrd2(const double *parTRAP, TString name)
1004 {
1005  // TRD2 cases
1006   Double_t *dummy=0;
1007   TString pbShape("BOX"), pbtiChonly("ONLY");
1008   if(name=="SCM0") {
1009     pbShape    = "TRD2";
1010     //    pbtiChonly = "MANY";
1011   }
1012   gMC->Gsvolu("PBTI", pbShape.Data(), fIdTmedArr[kIdPB], dummy, 0);
1013
1014   int nr=0;
1015   Info("PbInTrapForTrd2"," Pb tiles inside %s: shape %s :pbtiChonly %s\n nrstart %i\n", 
1016   name.Data(), pbShape.Data(), pbtiChonly.Data(), nr);
1017   AliEMCALGeometry * g = GetGeometry(); 
1018
1019   double par[5], parPB[5], parSC[5];
1020   //double fSampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
1021   double xpos = 0.0, ypos = 0.0;
1022   double zpos = -fSampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
1023   if(name == "SCMX") { // common trapezoid - 11 parameters
1024     double coef = (parTRAP[8] -  parTRAP[4]) / (2.*parTRAP[0]);
1025     double xCenterSCMX =  (parTRAP[4] +  parTRAP[8])/2.; // the same for y
1026     AliDebug(2,Form(" xCenterSCMX %8.5f : coef %8.7f \n", xCenterSCMX, coef));
1027
1028     par[2] = g->GetECPbRadThick()/2.; // z
1029     for(int iz=0; iz<g->GetNECLayers(); iz++){
1030       par[0] = parTRAP[4] + coef*fSampleWidth*iz;
1031       par[1] = par[0];
1032       xpos   = ypos = par[0] - xCenterSCMX;
1033     //if(parTRAP[1] < 0.) xpos = -xpos;
1034       gMC->Gsposp("PBTI", ++nr, name.Data(), xpos, ypos, zpos, 0, "ONLY", par, 3) ;
1035       AliDebug(2,Form(" %2.2i xpos %8.5f zpos %6.3f par[0,1] %6.3f |", iz+1, xpos, zpos, par[0]));
1036       if(iz%2>0) AliDebug(2,Form("\n"));
1037       zpos += fSampleWidth;
1038     } 
1039     AliDebug(2,Form(" Number of Pb tiles in SCMX %i coef %9.7f \n", nr, coef));
1040     AliDebug(2,Form(" par[1] %9.5f  par[2] %9.5f ypos %9.5f \n", par[1], par[2], ypos)); 
1041   } else if(name == "SCM0") { // 1-mar-05 ; TRD2 - 5 parameters
1042     AliDebug(2,Form(" SCM0 par = "));
1043     for(int i=0; i<5; i++) AliDebug(2,Form(" %9.5f ", parTRAP[i]));
1044     AliDebug(2,Form("\n zpos %f \n",zpos));
1045
1046     double tanx = (parTRAP[1] -  parTRAP[0]) / (2.*parTRAP[4]); //  tanx =  tany now
1047     double tany = (parTRAP[3] -  parTRAP[2]) / (2.*parTRAP[4]), ztmp=0.;
1048     parPB[4] = g->GetECPbRadThick()/2.;
1049     parSC[2] = g->GetECScintThick()/2.;
1050     for(int iz=0; iz<g->GetNECLayers(); iz++){
1051       ztmp     = fSampleWidth*double(iz);
1052       parPB[0] = parTRAP[0] + tanx*ztmp;
1053       parPB[1] = parPB[0]   + tanx*g->GetECPbRadThick();
1054       parPB[2] = parTRAP[2] + tany*ztmp;
1055       parPB[3] = parPB[2]   + tany*g->GetECPbRadThick();
1056       gMC->Gsposp("PBTI", ++nr, name.Data(), xpos, ypos, zpos, 0, pbtiChonly.Data(), parPB, 5) ;
1057       AliDebug(2,Form("\n PBTI %2i | zpos %6.3f | par = ", nr, zpos));
1058       /*
1059       for(int i=0; i<5; i++) printf(" %9.5f ", parPB[i]);
1060       // individual SC tile
1061       parSC[0] = parPB[0];
1062       parSC[1] = parPB[1];
1063       gMC->Gsposp("SCTI", nr, name.Data(), xpos, ypos, zpos+g->GetECScintThick(), 
1064       0, pbtiChonly.Data(), parSC, 3) ;
1065       printf("\n SCTI     zpos %6.3f | par = ", zpos+g->GetECScintThick());
1066       for(int i=0; i<3; i++) printf(" %9.5f ", parPB[i]);
1067       */
1068       zpos  += fSampleWidth;
1069     }
1070     AliDebug(2,Form("\n"));
1071   }
1072   Info("PbInTrapForTrd2", "Ver. 0.03 : was tested.");
1073 }
1074
1075 // 15-mar-05
1076 //______________________________________________________________________
1077 void AliEMCALv0::PbmoInTrd2(const AliEMCALGeometry * g, const Double_t emodPar[5], Double_t parPBMO[5])
1078 {
1079   // Pb inside Trd2
1080   Info("PbmoInTrd2"," started : geometry %s ", g->GetName());
1081   double wallThickness = g->GetPhiModuleSize()/2. -  g->GetPhiTileSize();
1082   AliDebug(2,Form(" wall thickness %7.5f \n", wallThickness));
1083   for(int i=0; i<4; i++) {
1084     parPBMO[i] = emodPar[i] - wallThickness;
1085     AliDebug(2,Form(" %i parPBMO %7.3f emodPar %7.3f : dif %7.3f \n", 
1086                     i, parPBMO[i],emodPar[i], parPBMO[i]-emodPar[i]));
1087   }
1088   parPBMO[4] = emodPar[4];
1089   gMC->Gsvolu("PBMO", "TRD2", fIdTmedArr[kIdPB], parPBMO, 5);
1090   gMC->Gspos("PBMO", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
1091   // Division 
1092   if(g->GetNPHIdiv()==2 && g->GetNETAdiv()==2) {
1093     Division2X2InPbmo(g, parPBMO);
1094     AliDebug(2,Form(" PBMO, division 2X2 | geometry |%s|\n", g->GetName()));
1095   } else {
1096     AliDebug(2,Form(" no division PBMO in this geometry |%s|\n", g->GetName()));
1097     assert(0);
1098   }
1099 }
1100
1101 //______________________________________________________________________
1102 void AliEMCALv0::Division2X2InPbmo(const AliEMCALGeometry * g, const Double_t parPBMO[5]) 
1103 {
1104   // Division 2X2
1105   Info("Division2X2InPbmo"," started : geometry %s ", g->GetName());
1106   //Double_t *dummy=0;
1107   //  gMC->Gsvolu("SCTI", "BOX", fIdTmedArr[kIdSC], dummy, 0);
1108
1109   double parSC[3];
1110   //double fSampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
1111   double xpos = 0.0, ypos = 0.0, zpos = 0.0, ztmp=0;;
1112   double tanx = (parPBMO[1] -  parPBMO[0]) / (2.*parPBMO[4]); //  tanx =  tany now
1113   double tany = (parPBMO[3] -  parPBMO[2]) / (2.*parPBMO[4]);
1114   char name[10], named[10], named2[10];
1115
1116   AliDebug(2,Form(" PBMO par = "));
1117   for(int i=0; i<5; i++) AliDebug(2,Form(" %9.5f ", parPBMO[i]));
1118   AliDebug(2,Form("\n"));
1119
1120   parSC[2] = g->GetECScintThick()/2.;
1121   zpos = -fSampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick() + g->GetECScintThick()/2.;
1122   AliDebug(2,Form(" parSC[2] %9.5f \n", parSC[2]));
1123   for(int iz=0; iz<g->GetNECLayers(); iz++){
1124     ztmp     = g->GetECPbRadThick() + fSampleWidth*double(iz); // Z for previous PB
1125     parSC[0] =  parPBMO[0] + tanx*ztmp;
1126     parSC[1] =  parPBMO[2] + tany*ztmp;
1127
1128     sprintf(name,"SC%2.2i", iz+1);
1129     gMC->Gsvolu(name, "BOX", fIdTmedArr[kIdSC], parSC, 3);
1130     gMC->Gspos(name, 1, "PBMO", xpos, ypos, zpos, 0, "ONLY") ;
1131     AliDebug(2,Form("%s | zpos %6.3f | parSC[0,1]=(%7.5f,%7.5f) -> ", 
1132                     name, zpos, parSC[0], parSC[1]));
1133     
1134     sprintf(named,"SY%2.2i", iz+1);
1135     printf(" %s -> ", named);
1136     gMC->Gsdvn(named,name, 2, 2);
1137
1138     sprintf(named2,"SX%2.2i", iz+1);
1139     printf(" %s \n", named2);
1140     gMC->Gsdvn(named2,named, 2, 1);
1141
1142     zpos    += fSampleWidth;
1143   }
1144 }
1145
1146 //______________________________________________________________________
1147 AliEMCALShishKebabTrd1Module* AliEMCALv0::GetShishKebabModule(Int_t neta)
1148
1149   // 28-oct-05
1150   AliEMCALShishKebabTrd1Module* trd1=0;
1151   if(fShishKebabModules && neta>=0 && neta<fShishKebabModules->GetSize()) {
1152     trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabModules->At(neta);
1153   }
1154   return trd1;
1155 }
1156
1157 //_____________________________________________________________________________
1158 void AliEMCALv0::AddAlignableVolumes() const
1159 {
1160   TString ntmp(GetTitle()); // name of EMCAL geometry
1161
1162   if(ntmp.Contains("WSUC")) {
1163     AddAlignableVolumesInWSUC();  // WSUC case
1164   } else {
1165     AddAlignableVolumesInALICE(); // ALICE case
1166   }
1167 }
1168
1169 //______________________________________________________________________
1170 void AliEMCALv0::AddAlignableVolumesInALICE() const
1171 {
1172   //
1173   // Create entries for alignable volumes associating the symbolic volume
1174   // name with the corresponding volume path. Needs to be synchronized with
1175   // eventual changes in the geometry.
1176   //
1177
1178   Float_t * pars = GetGeometry()->GetSuperModulesPars();
1179   double rpos = (GetGeometry()->GetEnvelop(0) + GetGeometry()->GetEnvelop(1))/2.;
1180   double phi, phiRad, xpos, ypos, zpos;
1181
1182   TString volpath, symname;
1183
1184   Int_t nSMod = GetGeometry()->GetNumberOfSuperModules(); 
1185   for (Int_t smodnum=0; smodnum < nSMod; smodnum++) {
1186     volpath = "ALIC_1/XEN1_1/SMOD_";
1187     volpath += (smodnum+1);
1188     symname = "EMCAL/FullSupermodule";
1189     symname += (smodnum+1);
1190
1191     if(GetGeometry()->GetKey110DEG() && smodnum>=10) {
1192       volpath = "ALIC_1/XEN1_1/SM10_";
1193       volpath += (smodnum-10+1);
1194       symname = "EMCAL/HalfSupermodule";    
1195       symname += (smodnum-10+1);
1196     }
1197
1198     if(!gGeoManager->SetAlignableEntry(symname.Data(),volpath.Data()))
1199       AliFatal("AliEMCALv0::Unable to set alignable entry!!");
1200
1201     // Creates the Tracking to Local transformation matrix for EMCAL
1202     // modules                         
1203     TGeoPNEntry *alignableEntry = gGeoManager->GetAlignableEntry(symname.Data()) ;
1204     const char *path = alignableEntry->GetTitle();
1205     if (!gGeoManager->cd(path))
1206       AliFatal(Form("Volume path %s not valid!",path));
1207
1208     phiRad = GetGeometry()->GetPhiCenterOfSM(smodnum);  //comes in radians, not degrees
1209     phi = phiRad*180./TMath::Pi();             //need degrees for rot. matrix
1210     xpos = rpos * TMath::Cos(phiRad);
1211     ypos = rpos * TMath::Sin(phiRad);
1212     zpos = pars[2];
1213     if(GetGeometry()->GetKey110DEG() && smodnum >= 10) {
1214       xpos += (pars[1]/2. * TMath::Sin(phiRad));
1215       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
1216     }
1217
1218     TGeoHMatrix *matTtoL;
1219     TGeoHMatrix* globMatrix = gGeoManager->GetCurrentMatrix();
1220
1221     if(smodnum%2 == 0) {
1222       // pozitive z                                                        
1223       TGeoTranslation geoTran0(xpos,ypos, zpos);                        
1224       TGeoRotation geoRot0("geoRot0", 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
1225       TGeoCombiTrans mat0(geoTran0, geoRot0);
1226       matTtoL = new TGeoHMatrix(mat0);
1227
1228       matTtoL->MultiplyLeft(&(globMatrix->Inverse()));
1229       alignableEntry->SetMatrix(matTtoL);
1230
1231     } else {
1232       // negative z                                                                                 
1233       double phiy = 90. + phi + 180.;
1234       if(phiy>=360.) phiy -= 360.;
1235       TGeoTranslation geoTran1(xpos,ypos,-zpos);
1236       TGeoRotation geoRot1("geoRot1", 90.0, phi, 90.0, phiy, 180.0, 0.0);
1237       TGeoCombiTrans mat1(geoTran1, geoRot1);
1238       matTtoL = new TGeoHMatrix(mat1);
1239
1240       matTtoL->MultiplyLeft(&(globMatrix->Inverse()));
1241       alignableEntry->SetMatrix(matTtoL);
1242
1243     }
1244
1245   }
1246
1247 }
1248
1249 //______________________________________________________________________
1250 void AliEMCALv0::AddAlignableVolumesInWSUC() const
1251 {
1252   //
1253   // Create entries for alignable volumes associating the symbolic volume
1254   // name with the corresponding volume path. Needs to be synchronized with
1255   // eventual changes in the geometry.
1256   //
1257
1258   TString vpstr1 = "WSUC_1/XEN1_1/SMOD_";
1259   TString snstr1 = "EMCAL/CosmicTestSupermodule";
1260   TString volpath, symname;
1261
1262   // #SM is just one 
1263   for (Int_t smodnum=0; smodnum < 1; smodnum++) {
1264     symname = snstr1;
1265     symname += (smodnum+1);
1266     volpath = vpstr1;
1267     volpath += (smodnum+1);
1268     if(!gGeoManager->SetAlignableEntry(symname.Data(),volpath.Data()))
1269       AliFatal("Unable to set alignable entry!!");
1270   }
1271 }