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, gAlice->GetMCApp()->GetCurrentTrackNumber(), iChamber, ipart,
 633                           fTrackPosition.X(), fTrackPosition.Y()+yAngleEffect, fTrackPosition.Z(), 0.0,
 634                           fTrackMomentum.P(),theta, phi, fStepSum[idvol], fDestepSum[idvol],
 635                           fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z());
 636 //     if (GetDebug()){
 637 //       Info("StepManager Exit","Particle exiting from chamber %d",iChamber);
 638 //       Info("StepManager Exit","StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]);
 639 //       Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 640 //     }
 641     fStepSum[idvol]  =0; // Reset for the next event
 642     fDestepSum[idvol]=0; // Reset for the next event
 643   }
 644 }
 645
 646 //__________________________________________
 647 void AliMUONv1::StepManagerOld()
 648 {
 649   // Old Stepmanager for the chambers
 650   Int_t          copy, id;
 651   static Int_t   idvol =-1;
 652   static Int_t   vol[2];
 653   Int_t          ipart;
 654   TLorentzVector pos;
 655   TLorentzVector mom;
 656   Float_t        theta,phi;
 657   Float_t        destep, step;
 658
 659   static Float_t sstep;
 660   static Float_t eloss, eloss2, xhit, yhit, zhit, tof, tlength;
 661   const  Float_t kBig = 1.e10;
 662   static Float_t hits[15];
 663
 664   TClonesArray &lhits = *fHits;
 665
 666   //
 667   //
 668   // Only charged tracks
 669   if( !(gMC->TrackCharge()) ) return;
 670   //
 671   // Only gas gap inside chamber
 672   // Tag chambers and record hits when track enters
 673   id=gMC->CurrentVolID(copy);
 674   vol[0] = GetChamberId(id);
 675   idvol = vol[0] -1;
 676
 677   if (idvol == -1) return;
 678
 679   //
 680   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 681   gMC->TrackPosition(pos);
 682   gMC->TrackMomentum(mom);
 683
 684   ipart  = gMC->TrackPid();
 685
 686   //
 687   // momentum loss and steplength in last step
 688   destep = gMC->Edep();
 689   step   = gMC->TrackStep();
 690   // cout<<"------------"<<step<<endl;
 691   //
 692   // record hits when track enters ...
 693   if( gMC->IsTrackEntering()) {
 694
 695       gMC->SetMaxStep(fMaxStepGas);
 696       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 697       Double_t rt = TMath::Sqrt(tc);
 698       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
 699       Double_t tx = mom[0]/pmom;
 700       Double_t ty = mom[1]/pmom;
 701       Double_t tz = mom[2]/pmom;
 702       Double_t s  = ((AliMUONChamber*)(*fChambers)[idvol])
 703           ->ResponseModel()
 704           ->Pitch()/tz;
 705       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 706       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 707       hits[0] = Float_t(ipart);         // Geant3 particle type
 708       hits[1] = pos[0]+s*tx;            // X-position for hit
 709       hits[2] = pos[1]+s*ty;            // Y-position for hit
 710       hits[3] = pos[2]+s*tz;            // Z-position for hit
 711       hits[4] = theta;                  // theta angle of incidence
 712       hits[5] = phi;                    // phi angle of incidence
 713       hits[8] = 0;//PadHits does not exist anymore  (Float_t) fNPadHits;    // first padhit
 714       hits[9] = -1;                     // last pad hit
 715       hits[10] = mom[3];                // hit momentum P
 716       hits[11] = mom[0];                // Px
 717       hits[12] = mom[1];                // Py
 718       hits[13] = mom[2];                // Pz
 719       tof=gMC->TrackTime();
 720       hits[14] = tof;                   // Time of flight
 721       tlength  = 0;
 722       eloss    = 0;
 723       eloss2   = 0;
 724       sstep=0;
 725       xhit     = pos[0];
 726       yhit     = pos[1];
 727       zhit     = pos[2];
 728       Chamber(idvol).ChargeCorrelationInit();
 729       // Only if not trigger chamber
 730
 731 //       printf("---------------------------\n");
 732 //       printf(">>>> Y =  %f \n",hits[2]);
 733 //       printf("---------------------------\n");
 734
 735
 736
 737      //  if(idvol < AliMUONConstants::NTrackingCh()) {
 738 //        //
 739 //        //  Initialize hit position (cursor) in the segmentation model
 740 //        ((AliMUONChamber*) (*fChambers)[idvol])
 741 //            ->SigGenInit(pos[0], pos[1], pos[2]);
 742 //       } else {
 743 //        //geant3->Gpcxyz();
 744 //        //printf("In the Trigger Chamber #%d\n",idvol-9);
 745 //       }
 746   }
 747   eloss2+=destep;
 748   sstep+=step;
 749
 750   // cout<<sstep<<endl;
 751
 752   //
 753   // Calculate the charge induced on a pad (disintegration) in case
 754   //
 755   // Mip left chamber ...
 756   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 757       gMC->SetMaxStep(kBig);
 758       eloss   += destep;
 759       tlength += step;
 760
 761       Float_t x0,y0,z0;
 762       Float_t localPos[3];
 763       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
 764       gMC->Gmtod(globalPos,localPos,1);
 765
 766       if(idvol < AliMUONConstants::NTrackingCh()) {
 767 // tracking chambers
 768           x0 = 0.5*(xhit+pos[0]);
 769           y0 = 0.5*(yhit+pos[1]);
 770           z0 = 0.5*(zhit+pos[2]);
 771       } else {
 772 // trigger chambers
 773           x0 = xhit;
 774           y0 = yhit;
 775           z0 = 0.;
 776       }
 777
 778
 779       //      if (eloss >0)  MakePadHits(x0,y0,z0,eloss,tof,idvol);
 780
 781
 782       hits[6] = tlength;   // track length
 783       hits[7] = eloss2;    // de/dx energy loss
 784
 785
 786       //      if (fNPadHits > (Int_t)hits[8]) {
 787       //          hits[8] = hits[8]+1;
 788       //          hits[9] = 0: // PadHits does not exist anymore (Float_t) fNPadHits;
 789       //}
 790 //
 791 //    new hit
 792
 793       new(lhits[fNhits++])
 794           AliMUONHit(fIshunt, gAlice->GetMCApp()->GetCurrentTrackNumber(), vol,hits);
 795       eloss = 0;
 796       //
 797       // Check additional signal generation conditions
 798       // defined by the segmentation
 799       // model (boundary crossing conditions)
 800       // only for tracking chambers
 801   } else if
 802       ((idvol < AliMUONConstants::NTrackingCh()) &&
 803        ((AliMUONChamber*) (*fChambers)[idvol])->SigGenCond(pos[0], pos[1], pos[2]))
 804   {
 805       ((AliMUONChamber*) (*fChambers)[idvol])
 806           ->SigGenInit(pos[0], pos[1], pos[2]);
 807
 808       Float_t localPos[3];
 809       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
 810       gMC->Gmtod(globalPos,localPos,1);
 811
 812       eloss    += destep;
 813
 814       // if (eloss > 0 && idvol < AliMUONConstants::NTrackingCh())
 815       //        MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),pos[2],eloss,tof,idvol);
 816       xhit     = pos[0];
 817       yhit     = pos[1];
 818       zhit     = pos[2];
 819       eloss = 0;
 820       tlength += step ;
 821       //
 822       // nothing special  happened, add up energy loss
 823   } else {
 824       eloss   += destep;
 825       tlength += step ;
 826   }
 827 }