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