Declaring flange as MANY.
[u/mrichter/AliRoot.git] / MUON / AliMUONv1.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *      SigmaEffect_thetadegrees                                                                  *
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 purpeateose. It is      *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 /////////////////////////////////////////////////////////
19 //  Manager and hits classes for set:MUON version 1    //
20 /////////////////////////////////////////////////////////
21
22 #include <TRandom.h>
23 #include <TF1.h>
24 #include <TClonesArray.h>
25 #include <TLorentzVector.h> 
26 #include <TNode.h> 
27 #include <TRandom.h> 
28 #include <TTUBE.h>
29 #include <TGeoMatrix.h>
30 #include <TVirtualMC.h>
31 #include <TParticle.h>
32
33 #include "AliConst.h" 
34 #include "AliMUONChamber.h"
35 #include "AliMUONConstants.h"
36 #include "AliMUONFactory.h"
37 #include "AliMUONHit.h"
38 #include "AliMUONTriggerCircuit.h"
39 #include "AliMUONv1.h"
40 #include "AliMUONVGeometryBuilder.h"    
41 #include "AliMUONChamberGeometry.h"     
42 #include "AliMUONGeometryEnvelope.h"    
43 #include "AliMUONGeometryConstituent.h" 
44 #include "AliMagF.h"
45 #include "AliRun.h"
46 #include "AliMC.h"
47
48 ClassImp(AliMUONv1)
49  
50 //___________________________________________
51 AliMUONv1::AliMUONv1() : AliMUON()
52   ,fTrackMomentum(), fTrackPosition(),fGlobalTransformation(0) 
53 {
54 // Constructor
55     fChambers   = 0;
56     fStepManagerVersionOld  = kFALSE;
57     fAngleEffect = kTRUE;
58     fStepMaxInActiveGas     = 0.6;
59     fStepSum    =  0x0;
60     fDestepSum  =  0x0;
61     fElossRatio =  0x0;
62     fAngleEffect10   = 0x0;
63     fAngleEffectNorma= 0x0;
64
65 //___________________________________________
66 AliMUONv1::AliMUONv1(const char *name, const char *title)
67   : AliMUON(name,title), fTrackMomentum(), fTrackPosition()
68 {
69 // Constructor
70     // By default include all stations
71     AliMUONFactory factory;
72     factory.Build(this, title);
73
74     fStepManagerVersionOld = kFALSE;
75     fAngleEffect = kTRUE;
76     fStepMaxInActiveGas = 0.6;
77
78     fStepSum   = new Float_t [AliMUONConstants::NCh()];
79     fDestepSum = new Float_t [AliMUONConstants::NCh()];
80     for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
81       fStepSum[i] =0.0;
82       fDestepSum[i]=0.0;
83     }
84     // Ratio of particle mean eloss with respect MIP's Khalil Boudjemline, sep 2003, PhD.Thesis and Particle Data Book
85     fElossRatio = new TF1("ElossRatio","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",0.5,5.); 
86     fElossRatio->SetParameter(0,1.02138);
87     fElossRatio->SetParameter(1,-9.54149e-02);
88     fElossRatio->SetParameter(2,+7.83433e-02); 
89     fElossRatio->SetParameter(3,-9.98208e-03);
90     fElossRatio->SetParameter(4,+3.83279e-04);
91
92     // Angle effect in tracking chambers at theta =10 degres as a function of ElossRatio (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) (in micrometers)
93     fAngleEffect10 = new TF1("AngleEffect10","[0]+[1]*x+[2]*x*x",0.5,3.0);
94     fAngleEffect10->SetParameter(0, 1.90691e+02);
95     fAngleEffect10->SetParameter(1,-6.62258e+01);
96     fAngleEffect10->SetParameter(2,+1.28247e+01);
97     // Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis)  
98     // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
99     fAngleEffectNorma = new TF1("AngleEffectNorma","[0]+[1]*x+[2]*x*x+[3]*x*x*x",0.0,10.0);
100     fAngleEffectNorma->SetParameter(0,4.148);
101     fAngleEffectNorma->SetParameter(1,-6.809e-01);
102     fAngleEffectNorma->SetParameter(2,5.151e-02);
103     fAngleEffectNorma->SetParameter(3,-1.490e-03);
104
105     // Define the global transformation:
106     // Transformation from the old ALICE coordinate system to a new one:
107     // x->-x, z->-z 
108     TGeoRotation* rotGlobal 
109       = new TGeoRotation("rotGlobal", 90., 180., 90., 90., 180., 0.);
110     fGlobalTransformation = new TGeoCombiTrans(0., 0., 0., rotGlobal);
111 }
112
113 //_____________________________________________________________________________
114 AliMUONv1::AliMUONv1(const AliMUONv1& right) 
115   : AliMUON(right) 
116 {  
117   // copy constructor (not implemented)
118
119   Fatal("AliMUONv1", "Copy constructor not provided.");
120 }
121
122 //___________________________________________
123 AliMUONv1::~AliMUONv1()
124 {
125 // Destructor
126
127   delete fGlobalTransformation;
128 }
129
130 //_____________________________________________________________________________
131 AliMUONv1& AliMUONv1::operator=(const AliMUONv1& right)
132 {
133   // assignement operator (not implemented)
134
135   // check assignement to self
136   if (this == &right) return *this;
137
138   Fatal("operator =", "Assignement operator not provided.");
139     
140   return *this;  
141 }    
142
143 //__________________________________________________
144 void AliMUONv1::CreateGeometry()
145 {
146 //
147 // Construct geometry using geometry builders.
148 //
149
150   for (Int_t i=0; i<fGeometryBuilders->GetEntriesFast(); i++) {
151
152     // Get the builder
153     AliMUONVGeometryBuilder* builder
154       = (AliMUONVGeometryBuilder*)fGeometryBuilders->At(i);
155
156     // Create geometry with each builder
157     if (builder) {
158       builder->CreateGeometry();
159       builder->SetTransformations();
160     }
161   }
162
163   for (Int_t j=0; j<AliMUONConstants::NCh(); j++) {
164
165     AliMUONChamberGeometry* geometry = Chamber(j).GetGeometry();
166
167     if (!geometry) continue;
168           // Skip chambers with not defined geometry  
169           
170     // Loop over envelopes
171     const TObjArray* kEnvelopes = geometry->GetEnvelopes();
172     for (Int_t k=0; k<kEnvelopes->GetEntriesFast(); k++) {
173
174       // Get envelope
175       AliMUONGeometryEnvelope* env = (AliMUONGeometryEnvelope*)kEnvelopes->At(k);
176       const TGeoCombiTrans* kEnvTrans = env->GetTransformation();
177       const char* only = "ONLY";
178       if (env->IsMANY()) only = "MANY";
179
180       if (env->IsVirtual() && env->GetConstituents()->GetEntriesFast() == 0 ) {
181         // virtual envelope + nof constituents = 0 
182         //         => not allowed;
183         //            empty virtual envelope has no sense 
184         Fatal("CreateGeometry", "Virtual envelope must have constituents.");
185         return;
186       }
187
188       if (!env->IsVirtual() && env->GetConstituents()->GetEntriesFast() > 0 ) {
189         // non virtual envelope + nof constituents > 0 
190         //        => not allowed;
191         //           use VMC to place constituents
192         Fatal("CreateGeometry", "Non virtual envelope cannot have constituents.");
193         return;
194       }
195
196       if (!env->IsVirtual() && env->GetConstituents()->GetEntriesFast() == 0 ) {
197         // non virtual envelope + nof constituents = 0 
198         //        => place envelope in ALICE by composed transformation:
199         //           Tglobal * Tch * Tenv
200
201         // Compound chamber transformation with the envelope one
202         TGeoHMatrix total 
203           = (*fGlobalTransformation) * 
204             (*geometry->GetTransformation()) * 
205             (*kEnvTrans);
206         PlaceVolume(env->GetName(), geometry->GetMotherVolume(),
207                     env->GetCopyNo(), total, 0, 0, only);
208       }
209
210       if (env->IsVirtual() && env->GetConstituents()->GetEntriesFast() > 0 ) {
211         // virtual envelope + nof constituents > 0 
212         //         => do not place envelope and place constituents
213         //            in ALICE by composed transformation:
214         //            Tglobal * Tch * Tenv * Tconst   
215
216         for  (Int_t l=0; l<env->GetConstituents()->GetEntriesFast(); l++) {
217           AliMUONGeometryConstituent* constituent
218             = (AliMUONGeometryConstituent*)env->GetConstituents()->At(l);
219
220           // Compound chamber transformation with the envelope one + the constituent one
221           TGeoHMatrix total 
222             = (*fGlobalTransformation) *
223               (*geometry->GetTransformation()) * 
224               (*kEnvTrans) * 
225               (*constituent->GetTransformation());
226
227           PlaceVolume(constituent->GetName(), geometry->GetMotherVolume(),
228                       constituent->GetCopyNo(), total,
229                       constituent->GetNpar(), constituent->GetParam(), only);
230         }
231       }
232     } 
233   }
234 }
235 //__________________________________________________________________
236 Int_t  AliMUONv1::GetChamberId(Int_t volId) const
237 {
238 // Check if the volume with specified  volId is a sensitive volume (gas) 
239 // of some chamber and returns the chamber number;
240 // if not sensitive volume - return 0.
241 // ---
242
243 /*
244   for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
245     if (volId==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) return i;
246 */
247   for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
248     if ( ((AliMUONChamber*)(*fChambers)[i-1])->IsSensId(volId) ) return i;
249
250   return 0;
251 }
252 //________________________________________________________________
253 void AliMUONv1::CreateMaterials()
254 {
255
256   // *** DEFINITION OF AVAILABLE MUON MATERIALS *** 
257   //
258   //     Ar-CO2 gas (80%+20%)
259     Float_t ag1[3]   = { 39.95,12.01,16. };
260     Float_t zg1[3]   = { 18.,6.,8. };
261     Float_t wg1[3]   = { .8,.0667,.13333 };
262     Float_t dg1      = .001821;
263     //
264     //     Ar-buthane-freon gas -- trigger chambers 
265     Float_t atr1[4]  = { 39.95,12.01,1.01,19. };
266     Float_t ztr1[4]  = { 18.,6.,1.,9. };
267     Float_t wtr1[4]  = { .56,.1262857,.2857143,.028 };
268     Float_t dtr1     = .002599;
269     //
270     //     Ar-CO2 gas 
271     Float_t agas[3]  = { 39.95,12.01,16. };
272     Float_t zgas[3]  = { 18.,6.,8. };
273     Float_t wgas[3]  = { .74,.086684,.173316 };
274     Float_t dgas     = .0018327;
275     //
276     //     Ar-Isobutane gas (80%+20%) -- tracking 
277     Float_t ag[3]    = { 39.95,12.01,1.01 };
278     Float_t zg[3]    = { 18.,6.,1. };
279     Float_t wg[3]    = { .8,.057,.143 };
280     Float_t dg       = .0019596;
281     //
282     //     Ar-Isobutane-Forane-SF6 gas (49%+7%+40%+4%) -- trigger 
283     Float_t atrig[5] = { 39.95,12.01,1.01,19.,32.066 };
284     Float_t ztrig[5] = { 18.,6.,1.,9.,16. };
285     Float_t wtrig[5] = { .49,1.08,1.5,1.84,0.04 };
286     Float_t dtrig    = .0031463;
287     //
288     //     bakelite 
289
290     Float_t abak[3] = {12.01 , 1.01 , 16.};
291     Float_t zbak[3] = {6.     , 1.   , 8.};
292     Float_t wbak[3] = {6.     , 6.   , 1.}; 
293     Float_t dbak = 1.4;
294
295     Float_t epsil, stmin, deemax, tmaxfd, stemax;
296
297     Int_t iSXFLD   = gAlice->Field()->Integ();
298     Float_t sXMGMX = gAlice->Field()->Max();
299     //
300     // --- Define the various materials for GEANT --- 
301     AliMaterial(9, "ALUMINIUM$", 26.98, 13., 2.7, 8.9, 37.2);
302     AliMaterial(10, "ALUMINIUM$", 26.98, 13., 2.7, 8.9, 37.2);
303     AliMaterial(15, "AIR$      ", 14.61, 7.3, .001205, 30423.24, 67500);
304     AliMixture(19, "Bakelite$", abak, zbak, dbak, -3, wbak);
305     AliMixture(20, "ArC4H10 GAS$", ag, zg, dg, 3, wg);
306     AliMixture(21, "TRIG GAS$", atrig, ztrig, dtrig, -5, wtrig);
307     AliMixture(22, "ArCO2 80%$", ag1, zg1, dg1, 3, wg1);
308     AliMixture(23, "Ar-freon $", atr1, ztr1, dtr1, 4, wtr1);
309     AliMixture(24, "ArCO2 GAS$", agas, zgas, dgas, 3, wgas);
310     // materials for slat: 
311     //     Sensitive area: gas (already defined) 
312     //     PCB: copper 
313     //     insulating material and frame: vetronite
314     //     walls: carbon, rohacell, carbon 
315   Float_t aglass[5]={12.01, 28.09, 16.,   10.8,  23.};
316   Float_t zglass[5]={ 6.,   14.,    8.,    5.,   11.};
317   Float_t wglass[5]={ 0.5,  0.105, 0.355, 0.03,  0.01};
318   Float_t dglass=1.74;
319
320   // rohacell: C9 H13 N1 O2
321   Float_t arohac[4] = {12.01,  1.01, 14.010, 16.};
322   Float_t zrohac[4] = { 6.,    1.,    7.,     8.};
323   Float_t wrohac[4] = { 9.,   13.,    1.,     2.};
324   Float_t drohac    = 0.03;
325
326   AliMaterial(31, "COPPER$",   63.54,    29.,   8.96,  1.4, 0.);
327   AliMixture(32, "Vetronite$",aglass, zglass, dglass,    5, wglass);
328   AliMaterial(33, "Carbon$",   12.01,     6.,  2.265, 18.8, 49.9);
329   AliMixture(34, "Rohacell$", arohac, zrohac, drohac,   -4, wrohac); 
330
331
332     epsil  = .001; // Tracking precision, 
333     stemax = -1.;  // Maximum displacement for multiple scat 
334     tmaxfd = -20.; // Maximum angle due to field deflection 
335     deemax = -.3;  // Maximum fractional energy loss, DLS 
336     stmin  = -.8;
337     //
338     //    Air 
339     AliMedium(1, "AIR_CH_US         ", 15, 1, iSXFLD, sXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
340     //
341     //    Aluminum 
342
343     AliMedium(4, "ALU_CH_US          ", 9, 0, iSXFLD, sXMGMX, tmaxfd, fMaxStepAlu, 
344             fMaxDestepAlu, epsil, stmin);
345     AliMedium(5, "ALU_CH_US          ", 10, 0, iSXFLD, sXMGMX, tmaxfd, fMaxStepAlu, 
346             fMaxDestepAlu, epsil, stmin);
347     //
348     //    Ar-isoC4H10 gas 
349
350     AliMedium(6, "AR_CH_US          ", 20, 1, iSXFLD, sXMGMX, tmaxfd, fMaxStepGas, 
351             fMaxDestepGas, epsil, stmin);
352 //
353     //    Ar-Isobuthane-Forane-SF6 gas 
354
355     AliMedium(7, "GAS_CH_TRIGGER    ", 21, 1, iSXFLD, sXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
356
357     AliMedium(8, "BAKE_CH_TRIGGER   ", 19, 0, iSXFLD, sXMGMX, tmaxfd, fMaxStepAlu, 
358             fMaxDestepAlu, epsil, stmin);
359
360     AliMedium(9, "ARG_CO2   ", 22, 1, iSXFLD, sXMGMX, tmaxfd, fMaxStepGas, 
361             fMaxDestepAlu, epsil, stmin);
362     // tracking media for slats: check the parameters!! 
363     AliMedium(11, "PCB_COPPER        ", 31, 0, iSXFLD, sXMGMX, tmaxfd, 
364               fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
365     AliMedium(12, "VETRONITE         ", 32, 0, iSXFLD, sXMGMX, tmaxfd, 
366               fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
367     AliMedium(13, "CARBON            ", 33, 0, iSXFLD, sXMGMX, tmaxfd, 
368               fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
369     AliMedium(14, "Rohacell          ", 34, 0, iSXFLD, sXMGMX, tmaxfd, 
370               fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
371
372
373
374   //.Materials specific to stations
375   // created via builders
376   
377   for (Int_t i=0; i<fGeometryBuilders->GetEntriesFast(); i++) {
378
379     // Get the builder
380     AliMUONVGeometryBuilder* builder
381       = (AliMUONVGeometryBuilder*)fGeometryBuilders->At(i);
382
383     // Create materials with each builder
384     if (builder) builder->CreateMaterials();
385   }
386 }
387
388 //______________________________________________________________________________
389 void AliMUONv1::PlaceVolume(const TString& name, const TString& mName, 
390                             Int_t copyNo, const TGeoHMatrix& matrix, 
391                             Int_t npar, Double_t* param, const char* only) const
392 {
393 // Place the volume specified by name with the given transformation matrix
394 // ---
395
396   // Do not apply global transformation 
397   // if mother volume == DDIP
398   // (as it is applied on this volume)
399   TGeoHMatrix transform(matrix);
400   if (mName == TString("DDIP")) {
401     transform = (*fGlobalTransformation) * transform;
402                // To be changed to (*fGlobalTransformation).inverse()
403                // when available in TGeo
404                // To make this correct also for a general case when
405                // (*fGlobalTransformation) * *fGlobalTransformation) != 1
406   }            
407      
408   // Decompose transformation
409   const Double_t* xyz = transform.GetTranslation();
410   const Double_t* rm = transform.GetRotationMatrix();
411         
412   //cout << "Got translation: "
413   //     << xyz[0] << " " << xyz[1] << " " << xyz[2] << endl;
414         
415   //cout << "Got rotation: "
416   //     << rm[0] << " " << rm[1] << " " << rm[2] << endl
417   //     << rm[3] << " " << rm[4] << " " << rm[5] << endl
418   //     << rm[6] << " " << rm[7] << " " << rm[8] << endl;
419
420   // Check for presence of rotation
421   // (will be nice to be available in TGeo)
422   const Double_t kTolerance = 1e-04;
423   Bool_t isRotation = true; 
424   if (TMath::Abs(rm[0] - 1.) < kTolerance &&
425       TMath::Abs(rm[1] - 0.) < kTolerance &&
426       TMath::Abs(rm[2] - 0.) < kTolerance &&
427       TMath::Abs(rm[3] - 0.) < kTolerance &&
428       TMath::Abs(rm[4] - 1.) < kTolerance &&
429       TMath::Abs(rm[5] - 0.) < kTolerance &&
430       TMath::Abs(rm[6] - 0.) < kTolerance &&
431       TMath::Abs(rm[7] - 0.) < kTolerance &&
432       TMath::Abs(rm[8] - 1.) < kTolerance) isRotation = false; 
433
434   Int_t krot = 0;
435   if (isRotation) {
436     TGeoRotation rot;
437     rot.SetMatrix(const_cast<Double_t*>(transform.GetRotationMatrix()));
438     Double_t theta1, phi1, theta2, phi2, theta3, phi3;
439     rot.GetAngles(theta1, phi1, theta2, phi2, theta3, phi3);
440         
441     //cout << "angles: " 
442     //     << theta1 << " " << phi1 << " "
443     //     << theta2 << " " << phi2 << " "
444     //     << theta3 << " " << phi3 << endl;
445         
446     AliMatrix(krot, theta1, phi1, theta2, phi2, theta3, phi3);
447   }     
448         
449   // Place the volume in ALIC
450   if (npar == 0)
451     gMC->Gspos(name, copyNo, mName, xyz[0], xyz[1], xyz[2] , krot, only);
452   else 
453     gMC->Gsposp(name, copyNo, mName, xyz[0], xyz[1], xyz[2] , krot, only,
454                 param, npar);
455
456
457
458 //___________________________________________
459 void AliMUONv1::Init()
460 {
461    // 
462    // Initialize Tracking Chambers
463    //
464
465    if(fDebug) printf("\n%s: Start Init for version 1 - CPC chamber type\n\n",ClassName());
466    Int_t i;
467    for (i=0; i<AliMUONConstants::NCh(); i++) {
468        ( (AliMUONChamber*) (*fChambers)[i])->Init();
469    }
470    
471    //
472    // Set the chamber (sensitive region) GEANT identifier
473    //
474    for (Int_t i=0; i<fGeometryBuilders->GetEntriesFast(); i++) {
475
476     // Get the builder
477     AliMUONVGeometryBuilder* builder
478       = (AliMUONVGeometryBuilder*)fGeometryBuilders->At(i);
479
480     // Set sesitive volumes with each builder
481     if (builder) builder->SetSensitiveVolumes();
482   }
483
484 /*
485    //
486    // Set the chamber (sensitive region) GEANT identifier
487    ((AliMUONChamber*)(*fChambers)[0])->SetGid(gMC->VolId("S01G"));
488    ((AliMUONChamber*)(*fChambers)[1])->SetGid(gMC->VolId("S02G"));
489
490    ((AliMUONChamber*)(*fChambers)[2])->SetGid(gMC->VolId("S03G"));
491    ((AliMUONChamber*)(*fChambers)[3])->SetGid(gMC->VolId("S04G"));
492
493    ((AliMUONChamber*)(*fChambers)[4])->SetGid(gMC->VolId("S05G"));
494    ((AliMUONChamber*)(*fChambers)[5])->SetGid(gMC->VolId("S06G"));
495
496    ((AliMUONChamber*)(*fChambers)[6])->SetGid(gMC->VolId("S07G"));
497    ((AliMUONChamber*)(*fChambers)[7])->SetGid(gMC->VolId("S08G"));
498
499    ((AliMUONChamber*)(*fChambers)[8])->SetGid(gMC->VolId("S09G"));
500    ((AliMUONChamber*)(*fChambers)[9])->SetGid(gMC->VolId("S10G"));
501
502    ((AliMUONChamber*)(*fChambers)[10])->SetGid(gMC->VolId("SG1A"));
503    ((AliMUONChamber*)(*fChambers)[11])->SetGid(gMC->VolId("SG2A"));
504    ((AliMUONChamber*)(*fChambers)[12])->SetGid(gMC->VolId("SG3A"));
505    ((AliMUONChamber*)(*fChambers)[13])->SetGid(gMC->VolId("SG4A"));
506 */
507    if(fDebug) printf("\n%s: Finished Init for version 1 - CPC chamber type\n",ClassName());
508
509    //cp 
510    if(fDebug) printf("\n%s: Start Init for Trigger Circuits\n",ClassName());
511    for (i=0; i<AliMUONConstants::NTriggerCircuit(); i++) {
512      ( (AliMUONTriggerCircuit*) (*fTriggerCircuits)[i])->Init(i);
513    }
514    if(fDebug) printf("%s: Finished Init for Trigger Circuits\n",ClassName());
515    //cp
516 }
517
518 //_______________________________________________________________________________
519 void AliMUONv1::StepManager()
520 {
521   // Stepmanager for the chambers
522
523  if (fStepManagerVersionOld) {
524     StepManagerOld();
525     return;
526   }
527
528   // Only charged tracks
529   if( !(gMC->TrackCharge()) ) return; 
530   // Only charged tracks
531   
532   // Only gas gap inside chamber
533   // Tag chambers and record hits when track enters 
534   static Int_t   idvol=-1;
535   Int_t   iChamber=0;
536   Int_t   id=0;
537   Int_t   copy;
538   const  Float_t kBig = 1.e10;
539
540
541   //
542   // Only gas gap inside chamber
543   // Tag chambers and record hits when track enters 
544   id=gMC->CurrentVolID(copy);
545   iChamber = GetChamberId(id);
546   idvol = iChamber -1;
547
548   if (idvol == -1) return;
549
550    if( gMC->IsTrackEntering() ) {
551      Float_t theta = fTrackMomentum.Theta();
552      if ((TMath::Pi()-theta)*kRaddeg>=15.) gMC->SetMaxStep(fStepMaxInActiveGas); // We use Pi-theta because z is negative
553   }
554
555 //  if (GetDebug()) {
556 //     Float_t z = ( (AliMUONChamber*)(*fChambers)[idvol])->Z() ;
557 //      Info("StepManager Step","Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber, z);
558 //   }  
559   // Particule id and mass, 
560   Int_t     ipart = gMC->TrackPid();
561   Float_t   mass  = gMC->TrackMass();
562
563   fDestepSum[idvol]+=gMC->Edep();
564   // Get current particle id (ipart), track position (pos)  and momentum (mom)
565   if ( fStepSum[idvol]==0.0 )  gMC->TrackMomentum(fTrackMomentum);
566   fStepSum[idvol]+=gMC->TrackStep();
567   
568 //   if (GetDebug()) {
569 //     Info("StepManager Step","iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g",
570 //       iChamber,ipart, fTrackMomentum.Theta()*kRaddeg, fTrackMomentum.Phi()*kRaddeg, mass, fStepSum[idvol], gMC->Edep());
571 //     Info("StepManager Step","Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(), fTrackMomentum.Z()) ;
572 //     gMC->TrackPosition(fTrackPosition);
573 //     Info("StepManager Step","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
574 //   }
575
576   // Track left chamber or StepSum larger than fStepMaxInActiveGas
577   if ( gMC->IsTrackExiting() || 
578        gMC->IsTrackStop() || 
579        gMC->IsTrackDisappeared()||
580        (fStepSum[idvol]>fStepMaxInActiveGas) ) {
581     
582     if   ( gMC->IsTrackExiting() || 
583            gMC->IsTrackStop() || 
584            gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
585
586     gMC->TrackPosition(fTrackPosition);
587     Float_t theta = fTrackMomentum.Theta();
588     Float_t phi   = fTrackMomentum.Phi();
589     
590     TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
591                                fStepSum[idvol]/2.*sin(theta)*sin(phi),
592                                fStepSum[idvol]/2.*cos(theta),0.0       );
593     //     if (GetDebug()) 
594     //       Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
595     //     if (GetDebug()) 
596     //        Info("StepManager Exit ","Track backToWire %f %f %f",backToWire.X(),backToWire.Y(),backToWire.Z()) ;
597     fTrackPosition-=backToWire;
598     
599     //-------------- Angle effect 
600     // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
601     
602     Float_t BetaxGamma    = fTrackMomentum.P()/mass;//  pc/mc2
603     Float_t sigmaEffect10degrees;
604     Float_t sigmaEffectThetadegrees;
605     Float_t eLossParticleELossMip;
606     Float_t yAngleEffect=0.;
607     Float_t thetawires      =  TMath::Abs( TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) );// We use Pi-theta because z is negative
608
609
610     if (fAngleEffect){
611     if ( (BetaxGamma >3.2)   &&  (thetawires*kRaddeg<=15.) ) {
612       BetaxGamma=TMath::Log(BetaxGamma);
613       eLossParticleELossMip = fElossRatio->Eval(BetaxGamma);
614       // 10 degrees is a reference for a model (arbitrary)
615       sigmaEffect10degrees=fAngleEffect10->Eval(eLossParticleELossMip);// in micrometers
616       // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
617       sigmaEffectThetadegrees =  sigmaEffect10degrees/fAngleEffectNorma->Eval(thetawires*kRaddeg);  // For 5mm gap  
618       if ( (iChamber==1)  ||  (iChamber==2) )  
619         sigmaEffectThetadegrees/=(1.09833e+00+1.70000e-02*(thetawires*kRaddeg)); // The gap is different (4mm)
620       yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
621     }
622     }
623     
624     // One hit per chamber
625     GetMUONData()->AddHit(fIshunt, gAlice->GetMCApp()->GetCurrentTrackNumber(), iChamber, ipart, 
626                           fTrackPosition.X(), fTrackPosition.Y()+yAngleEffect, fTrackPosition.Z(), 0.0, 
627                           fTrackMomentum.P(),theta, phi, fStepSum[idvol], fDestepSum[idvol],
628                           fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z());
629 //     if (GetDebug()){
630 //       Info("StepManager Exit","Particle exiting from chamber %d",iChamber);
631 //       Info("StepManager Exit","StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]);
632 //       Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
633 //     }
634     fStepSum[idvol]  =0; // Reset for the next event
635     fDestepSum[idvol]=0; // Reset for the next event
636   }
637 }
638
639 //__________________________________________
640 void AliMUONv1::StepManagerOld()
641 {
642   // Old Stepmanager for the chambers
643   Int_t          copy, id;
644   static Int_t   idvol =-1;
645   static Int_t   vol[2];
646   Int_t          ipart;
647   TLorentzVector pos;
648   TLorentzVector mom;
649   Float_t        theta,phi;
650   Float_t        destep, step;
651   
652   static Float_t sstep;
653   static Float_t eloss, eloss2, xhit, yhit, zhit, tof, tlength;
654   const  Float_t kBig = 1.e10;
655   static Float_t hits[15];
656
657   TClonesArray &lhits = *fHits;
658
659   //
660   //
661   // Only charged tracks
662   if( !(gMC->TrackCharge()) ) return; 
663   //
664   // Only gas gap inside chamber
665   // Tag chambers and record hits when track enters 
666   id=gMC->CurrentVolID(copy);
667   vol[0] = GetChamberId(id);
668   idvol = vol[0] -1;
669
670   if (idvol == -1) return;
671
672   //
673   // Get current particle id (ipart), track position (pos)  and momentum (mom) 
674   gMC->TrackPosition(pos);
675   gMC->TrackMomentum(mom);
676
677   ipart  = gMC->TrackPid();
678
679   //
680   // momentum loss and steplength in last step
681   destep = gMC->Edep();
682   step   = gMC->TrackStep();
683   // cout<<"------------"<<step<<endl;
684   //
685   // record hits when track enters ...
686   if( gMC->IsTrackEntering()) {
687
688       gMC->SetMaxStep(fMaxStepGas);
689       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
690       Double_t rt = TMath::Sqrt(tc);
691       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
692       Double_t tx = mom[0]/pmom;
693       Double_t ty = mom[1]/pmom;
694       Double_t tz = mom[2]/pmom;
695       Double_t s  = ((AliMUONChamber*)(*fChambers)[idvol])
696           ->ResponseModel()
697           ->Pitch()/tz;
698       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
699       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
700       hits[0] = Float_t(ipart);         // Geant3 particle type
701       hits[1] = pos[0]+s*tx;            // X-position for hit
702       hits[2] = pos[1]+s*ty;            // Y-position for hit
703       hits[3] = pos[2]+s*tz;            // Z-position for hit
704       hits[4] = theta;                  // theta angle of incidence
705       hits[5] = phi;                    // phi angle of incidence 
706       hits[8] = 0;//PadHits does not exist anymore  (Float_t) fNPadHits;    // first padhit
707       hits[9] = -1;                     // last pad hit
708       hits[10] = mom[3];                // hit momentum P
709       hits[11] = mom[0];                // Px
710       hits[12] = mom[1];                // Py
711       hits[13] = mom[2];                // Pz
712       tof=gMC->TrackTime();
713       hits[14] = tof;                   // Time of flight
714       tlength  = 0;
715       eloss    = 0;
716       eloss2   = 0;
717       sstep=0;
718       xhit     = pos[0];
719       yhit     = pos[1];      
720       zhit     = pos[2];      
721       Chamber(idvol).ChargeCorrelationInit();
722       // Only if not trigger chamber
723
724 //       printf("---------------------------\n");
725 //       printf(">>>> Y =  %f \n",hits[2]);
726 //       printf("---------------------------\n");
727     
728       
729
730      //  if(idvol < AliMUONConstants::NTrackingCh()) {
731 //        //
732 //        //  Initialize hit position (cursor) in the segmentation model 
733 //        ((AliMUONChamber*) (*fChambers)[idvol])
734 //            ->SigGenInit(pos[0], pos[1], pos[2]);
735 //       } else {
736 //        //geant3->Gpcxyz();
737 //        //printf("In the Trigger Chamber #%d\n",idvol-9);
738 //       }
739   }
740   eloss2+=destep;
741   sstep+=step;
742
743   // cout<<sstep<<endl;
744
745   // 
746   // Calculate the charge induced on a pad (disintegration) in case 
747   //
748   // Mip left chamber ...
749   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
750       gMC->SetMaxStep(kBig);
751       eloss   += destep;
752       tlength += step;
753       
754       Float_t x0,y0,z0;
755       Float_t localPos[3];
756       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
757       gMC->Gmtod(globalPos,localPos,1); 
758
759       if(idvol < AliMUONConstants::NTrackingCh()) {
760 // tracking chambers
761           x0 = 0.5*(xhit+pos[0]);
762           y0 = 0.5*(yhit+pos[1]);
763           z0 = 0.5*(zhit+pos[2]);
764       } else {
765 // trigger chambers
766           x0 = xhit;
767           y0 = yhit;
768           z0 = 0.;
769       }
770       
771
772       //      if (eloss >0)  MakePadHits(x0,y0,z0,eloss,tof,idvol);
773       
774           
775       hits[6] = tlength;   // track length
776       hits[7] = eloss2;    // de/dx energy loss
777
778
779       //      if (fNPadHits > (Int_t)hits[8]) {
780       //          hits[8] = hits[8]+1;
781       //          hits[9] = 0: // PadHits does not exist anymore (Float_t) fNPadHits;
782       //}
783 //
784 //    new hit 
785       
786       new(lhits[fNhits++]) 
787           AliMUONHit(fIshunt, gAlice->GetMCApp()->GetCurrentTrackNumber(), vol,hits);
788       eloss = 0; 
789       //
790       // Check additional signal generation conditions 
791       // defined by the segmentation
792       // model (boundary crossing conditions)
793       // only for tracking chambers
794   } else if 
795       ((idvol < AliMUONConstants::NTrackingCh()) &&
796        ((AliMUONChamber*) (*fChambers)[idvol])->SigGenCond(pos[0], pos[1], pos[2]))
797   {
798       ((AliMUONChamber*) (*fChambers)[idvol])
799           ->SigGenInit(pos[0], pos[1], pos[2]);
800       
801       Float_t localPos[3];
802       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
803       gMC->Gmtod(globalPos,localPos,1); 
804
805       eloss    += destep;
806
807       // if (eloss > 0 && idvol < AliMUONConstants::NTrackingCh())
808       //        MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),pos[2],eloss,tof,idvol);
809       xhit     = pos[0];
810       yhit     = pos[1]; 
811       zhit     = pos[2];
812       eloss = 0;
813       tlength += step ;
814       //
815       // nothing special  happened, add up energy loss
816   } else {        
817       eloss   += destep;
818       tlength += step ;
819   }
820 }