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