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 <TVirtualMC.h>
  27 #include <TGeoMatrix.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 "AliMUONGeometryBuilder.h"
  37 #include "AliMUONGeometryModule.h"
  38 #include "AliMUONGeometrySVMap.h"
  39 #include "AliMUONGeometryDetElement.h"
  40 #include "AliMagF.h"
  41 #include "AliRun.h"
  42 #include "AliMC.h"
  43 #include "AliLog.h"
  44
  45 ClassImp(AliMUONv1)
  46
  47 //___________________________________________
  48 AliMUONv1::AliMUONv1()
  49   : AliMUON(),
  50     fStepManagerVersionOld(kFALSE),
  51     fStepManagerVersionDE(kFALSE),
  52     fAngleEffect(kTRUE),
  53     fStepMaxInActiveGas(0.6),
  54     fStepSum(0x0),
  55     fDestepSum(0x0),
  56     fTrackMomentum(),
  57     fTrackPosition(),
  58     fElossRatio(0x0),
  59     fAngleEffect10(0x0),
  60     fAngleEffectNorma(0x0)
  61 {
  62 // Default constructor
  63 }
  64
  65 //___________________________________________
  66 AliMUONv1::AliMUONv1(const char *name, const char *title)
  67   : AliMUON(name,title),
  68     fStepManagerVersionOld(kFALSE),
  69     fStepManagerVersionDE(kFALSE),
  70     fAngleEffect(kTRUE),
  71     fStepMaxInActiveGas(0.6),
  72     fStepSum(0x0),
  73     fDestepSum(0x0),
  74     fTrackMomentum(),
  75     fTrackPosition(),
  76     fElossRatio(0x0),
  77     fAngleEffect10(0x0),
  78     fAngleEffectNorma(0x0)
  79 {
  80 // Standard onstructor
  81
  82     // By default include all stations
  83     AliMUONFactory factory;
  84     factory.Build(this, title);
  85
  86     fStepSum   = new Float_t [AliMUONConstants::NCh()];
  87     fDestepSum = new Float_t [AliMUONConstants::NCh()];
  88     for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
  89       fStepSum[i] =0.0;
  90       fDestepSum[i]=0.0;
  91     }
  92     // Ratio of particle mean eloss with respect MIP's Khalil Boudjemline, sep 2003, PhD.Thesis and Particle Data Book
  93     fElossRatio = new TF1("ElossRatio","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",0.5,5.);
  94     fElossRatio->SetParameter(0,1.02138);
  95     fElossRatio->SetParameter(1,-9.54149e-02);
  96     fElossRatio->SetParameter(2,+7.83433e-02);
  97     fElossRatio->SetParameter(3,-9.98208e-03);
  98     fElossRatio->SetParameter(4,+3.83279e-04);
  99
 100     // Angle effect in tracking chambers at theta =10 degres as a function of ElossRatio (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) (in micrometers)
 101     fAngleEffect10 = new TF1("AngleEffect10","[0]+[1]*x+[2]*x*x",0.5,3.0);
 102     fAngleEffect10->SetParameter(0, 1.90691e+02);
 103     fAngleEffect10->SetParameter(1,-6.62258e+01);
 104     fAngleEffect10->SetParameter(2,+1.28247e+01);
 105     // Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis)
 106     // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
 107     fAngleEffectNorma = new TF1("AngleEffectNorma","[0]+[1]*x+[2]*x*x+[3]*x*x*x",0.0,10.0);
 108     fAngleEffectNorma->SetParameter(0,4.148);
 109     fAngleEffectNorma->SetParameter(1,-6.809e-01);
 110     fAngleEffectNorma->SetParameter(2,5.151e-02);
 111     fAngleEffectNorma->SetParameter(3,-1.490e-03);
 112 }
 113
 114 //_____________________________________________________________________________
 115 AliMUONv1::AliMUONv1(const AliMUONv1& right)
 116   : AliMUON(right)
 117 {
 118   // copy constructor (not implemented)
 119
 120   AliFatal("Copy constructor not provided.");
 121 }
 122
 123 //___________________________________________
 124 AliMUONv1::~AliMUONv1()
 125 {
 126 // Destructor
 127
 128   delete [] fStepSum;
 129   delete [] fDestepSum;
 130   delete fElossRatio;
 131   delete fAngleEffect10;
 132   delete fAngleEffectNorma;
 133 }
 134
 135 //_____________________________________________________________________________
 136 AliMUONv1& AliMUONv1::operator=(const AliMUONv1& right)
 137 {
 138   // assignement operator (not implemented)
 139
 140   // check assignement to self
 141   if (this == &right) return *this;
 142
 143   AliFatal("Assignement operator not provided.");
 144
 145   return *this;
 146 }
 147
 148 //__________________________________________________
 149 void AliMUONv1::CreateGeometry()
 150 {
 151 //
 152 // Construct geometry using geometry builder
 153 //
 154
 155   fGeometryBuilder->CreateGeometry();
 156 }
 157
 158 //________________________________________________________________
 159 void AliMUONv1::CreateMaterials()
 160 {
 161 //
 162 // Construct materials using geometry builder
 163 //
 164
 165   fGeometryBuilder->CreateMaterials();
 166 }
 167
 168 //___________________________________________
 169 void AliMUONv1::Init()
 170 {
 171    //
 172    // Initialize Tracking Chambers
 173    //
 174    AliDebug(1,"Start Init for version 1 - CPC chamber type");
 175    Int_t i;
 176    for (i=0; i<AliMUONConstants::NCh(); i++) {
 177        ( (AliMUONChamber*) (*fChambers)[i])->Init();
 178    }
 179
 180    //
 181    // Initialize geometry
 182    //
 183    fGeometryBuilder->InitGeometry();
 184    AliDebug(1,"Finished Init for version 1 - CPC chamber type");
 185
 186    //cp
 187    AliDebug(1,"Start Init for Trigger Circuits");
 188    for (i=0; i<AliMUONConstants::NTriggerCircuit(); i++) {
 189      ( (AliMUONTriggerCircuit*) (*fTriggerCircuits)[i])->Init(i);
 190    }
 191    AliDebug(1,"Finished Init for Trigger Circuits");
 192    //cp
 193 }
 194
 195 //__________________________________________________________________
 196 Int_t  AliMUONv1::GetChamberId(Int_t volId) const
 197 {
 198 // Check if the volume with specified  volId is a sensitive volume (gas)
 199 // of some chamber and returns the chamber number;
 200 // if not sensitive volume - return 0.
 201 // ---
 202
 203 /*
 204   for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
 205     if (volId==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) return i;
 206 */
 207   for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
 208     if ( ((AliMUONChamber*)(*fChambers)[i-1])->IsSensId(volId) ) return i;
 209
 210   return 0;
 211 }
 212
 213 //_______________________________________________________________________________
 214 TString  AliMUONv1::CurrentVolumePath() const
 215 {
 216 // Returns current volume path
 217 // (Could be removed when this function is available via gMC)
 218 // ---
 219
 220   TString path = "";
 221   TString name;
 222   Int_t copyNo;
 223   Int_t imother = 0;
 224   do {
 225     name = gMC->CurrentVolOffName(imother);
 226     gMC->CurrentVolOffID(imother++, copyNo);
 227     TString add = "/";
 228     add += name;
 229     add += ".";
 230     add += copyNo;
 231     path.Insert(0,add);
 232   }
 233   while ( name != TString("ALIC") );
 234
 235   return path;
 236 }
 237
 238 //_______________________________________________________________________________
 239 void AliMUONv1::StepManager()
 240 {
 241   // Stepmanager for the chambers
 242   // TBR
 243
 244  if (fStepManagerVersionDE) {
 245     StepManager2();
 246     return;
 247   }
 248
 249  if (fStepManagerVersionOld) {
 250     StepManagerOld();
 251     return;
 252   }
 253
 254   // Only charged tracks
 255   if( !(gMC->TrackCharge()) ) return;
 256   // Only charged tracks
 257
 258   // Only gas gap inside chamber
 259   // Tag chambers and record hits when track enters
 260   static Int_t   idvol=-1;
 261   Int_t   iChamber=0;
 262   Int_t   id=0;
 263   Int_t   copy;
 264   const  Float_t kBig = 1.e10;
 265
 266
 267   //
 268   // Only gas gap inside chamber
 269   // Tag chambers and record hits when track enters
 270   id=gMC->CurrentVolID(copy);
 271   iChamber = GetChamberId(id);
 272   idvol = iChamber -1;
 273
 274   if (idvol == -1) return;
 275
 276   // Filling TrackRefs file for MUON. Our Track references are the active volume of the chambers
 277   if ( (gMC->IsTrackEntering() || gMC->IsTrackExiting() ) )
 278     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 279
 280    if( gMC->IsTrackEntering() ) {
 281      Float_t theta = fTrackMomentum.Theta();
 282      if ((TMath::Pi()-theta)*kRaddeg>=15.) gMC->SetMaxStep(fStepMaxInActiveGas); // We use Pi-theta because z is negative
 283   }
 284
 285 //  if (GetDebug()) {
 286 //     Float_t z = ( (AliMUONChamber*)(*fChambers)[idvol])->Z() ;
 287 //      Info("StepManager Step","Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber, z);
 288 //   }
 289   // Particule id and mass,
 290   Int_t     ipart = gMC->TrackPid();
 291   Float_t   mass  = gMC->TrackMass();
 292
 293   fDestepSum[idvol]+=gMC->Edep();
 294   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 295   if ( fStepSum[idvol]==0.0 )  gMC->TrackMomentum(fTrackMomentum);
 296   fStepSum[idvol]+=gMC->TrackStep();
 297
 298 //   if (GetDebug()) {
 299 //     Info("StepManager Step","iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g",
 300 //       iChamber,ipart, fTrackMomentum.Theta()*kRaddeg, fTrackMomentum.Phi()*kRaddeg, mass, fStepSum[idvol], gMC->Edep());//     Info("StepManager Step","Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(), fTrackMomentum.Z()) ;
 301 //     gMC->TrackPosition(fTrackPosition);
 302 //     Info("StepManager Step","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 303 //   }
 304
 305   // Track left chamber or StepSum larger than fStepMaxInActiveGas
 306   if ( gMC->IsTrackExiting() ||
 307        gMC->IsTrackStop() ||
 308        gMC->IsTrackDisappeared()||
 309        (fStepSum[idvol]>fStepMaxInActiveGas) ) {
 310
 311     if   ( gMC->IsTrackExiting() ||
 312            gMC->IsTrackStop() ||
 313            gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
 314
 315     gMC->TrackPosition(fTrackPosition);
 316     Float_t theta = fTrackMomentum.Theta();
 317     Float_t phi   = fTrackMomentum.Phi();
 318
 319     TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
 320                                fStepSum[idvol]/2.*sin(theta)*sin(phi),
 321                                fStepSum[idvol]/2.*cos(theta),0.0       );
 322     //     if (GetDebug())
 323     //       Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 324     //     if (GetDebug())
 325     //        Info("StepManager Exit ","Track backToWire %f %f %f",backToWire.X(),backToWire.Y(),backToWire.Z()) ;
 326     fTrackPosition-=backToWire;
 327
 328     //-------------- Angle effect
 329     // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
 330
 331     Float_t betaxGamma    = fTrackMomentum.P()/mass;//  pc/mc2
 332     Float_t sigmaEffect10degrees;
 333     Float_t sigmaEffectThetadegrees;
 334     Float_t eLossParticleELossMip;
 335     Float_t yAngleEffect=0.;
 336     Float_t thetawires      =  TMath::Abs( TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) );// We use Pi-theta because z is negative
 337
 338
 339     if (fAngleEffect){
 340     if ( (betaxGamma >3.2)   &&  (thetawires*kRaddeg<=15.) ) {
 341       betaxGamma=TMath::Log(betaxGamma);
 342       eLossParticleELossMip = fElossRatio->Eval(betaxGamma);
 343       // 10 degrees is a reference for a model (arbitrary)
 344       sigmaEffect10degrees=fAngleEffect10->Eval(eLossParticleELossMip);// in micrometers
 345       // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
 346       sigmaEffectThetadegrees =  sigmaEffect10degrees/fAngleEffectNorma->Eval(thetawires*kRaddeg);  // For 5mm gap
 347       if ( (iChamber==1)  ||  (iChamber==2) )
 348         sigmaEffectThetadegrees/=(1.09833e+00+1.70000e-02*(thetawires*kRaddeg)); // The gap is different (4mm)
 349       yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
 350     }
 351     }
 352
 353     // Detection elements ids
 354     AliMUONGeometryModule* geometry
 355       = Chamber(iChamber-1).GetGeometry();
 356
 357     AliMUONGeometryDetElement* detElement
 358       = geometry->FindBySensitiveVolume(CurrentVolumePath());
 359
 360     Int_t detElemId = 0;
 361     if (detElement) detElemId = detElement->GetUniqueID();
 362
 363     if (!detElemId) {
 364       cerr << "Chamber id: "
 365            << setw(3) << iChamber << "  "
 366            << "Current SV: "
 367            <<  CurrentVolumePath()
 368            << "  detElemId: "
 369            << setw(5) << detElemId
 370            << endl;
 371       Double_t x, y, z;
 372       gMC->TrackPosition(x, y, z);
 373       cerr << "   global position: "
 374            << x << ", " << y << ", " << z
 375            << endl;
 376       AliWarning("DetElemId not identified.");
 377     }
 378
 379     // One hit per chamber
 380     GetMUONData()->AddHit(fIshunt,
 381                           gAlice->GetMCApp()->GetCurrentTrackNumber(),
 382                           iChamber, ipart,
 383                           fTrackPosition.X(),
 384                           fTrackPosition.Y()+yAngleEffect,
 385                           fTrackPosition.Z(),
 386                           gMC->TrackTime(),
 387                           fTrackMomentum.P(),
 388                           theta,
 389                           phi,
 390                           fStepSum[idvol],
 391                           fDestepSum[idvol],
 392                           fTrackPosition.X(),
 393                           fTrackPosition.Y(),
 394                           fTrackPosition.Z());
 395
 396 //     if (GetDebug()){
 397 //       Info("StepManager Exit","Particle exiting from chamber %d",iChamber);
 398 //       Info("StepManager Exit","StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]);
 399 //       Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 400 //     }
 401     fStepSum[idvol]  =0; // Reset for the next event
 402     fDestepSum[idvol]=0; // Reset for the next event
 403   }
 404 }
 405
 406 //_______________________________________________________________________________
 407 void AliMUONv1::StepManager2()
 408 {
 409   // Stepmanager for the chambers
 410
 411  if (fStepManagerVersionOld) {
 412     StepManagerOld2();
 413     return;
 414   }
 415
 416   // Only charged tracks
 417   if( !(gMC->TrackCharge()) ) return;
 418   // Only charged tracks
 419
 420   // Only gas gap inside chamber
 421   // Tag chambers and record hits when track enters
 422   static Int_t   idvol=-1;
 423   Int_t   iChamber=0;
 424   Int_t   id=0;
 425   Int_t   copy;
 426   const  Float_t kBig = 1.e10;
 427
 428
 429   //
 430   // Only gas gap inside chamber
 431   // Tag chambers and record hits when track enters
 432   id=gMC->CurrentVolID(copy);
 433   iChamber = GetChamberId(id);
 434   idvol = iChamber -1;
 435
 436   if (idvol == -1) return;
 437
 438   // Filling TrackRefs file for MUON. Our Track references are the active volume of the chambers
 439   if ( (gMC->IsTrackEntering() || gMC->IsTrackExiting() ) )
 440     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 441
 442    if( gMC->IsTrackEntering() ) {
 443      Float_t theta = fTrackMomentum.Theta();
 444      if ((TMath::Pi()-theta)*kRaddeg>=15.) gMC->SetMaxStep(fStepMaxInActiveGas); // We use Pi-theta because z is negative
 445   }
 446
 447 //  if (GetDebug()) {
 448 //     Float_t z = ( (AliMUONChamber*)(*fChambers)[idvol])->Z() ;
 449 //      Info("StepManager Step","Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber, z);
 450 //   }
 451   // Particule id and mass,
 452   Int_t     ipart = gMC->TrackPid();
 453   Float_t   mass  = gMC->TrackMass();
 454
 455   fDestepSum[idvol]+=gMC->Edep();
 456   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 457   if ( fStepSum[idvol]==0.0 )  gMC->TrackMomentum(fTrackMomentum);
 458   fStepSum[idvol]+=gMC->TrackStep();
 459
 460 //   if (GetDebug()) {
 461 //     Info("StepManager Step","iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g",
 462 //       iChamber,ipart, fTrackMomentum.Theta()*kRaddeg, fTrackMomentum.Phi()*kRaddeg, mass, fStepSum[idvol], gMC->Edep());//     Info("StepManager Step","Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(), fTrackMomentum.Z()) ;
 463 //     gMC->TrackPosition(fTrackPosition);
 464 //     Info("StepManager Step","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 465 //   }
 466
 467   // Track left chamber or StepSum larger than fStepMaxInActiveGas
 468   if ( gMC->IsTrackExiting() ||
 469        gMC->IsTrackStop() ||
 470        gMC->IsTrackDisappeared()||
 471        (fStepSum[idvol]>fStepMaxInActiveGas) ) {
 472
 473     if   ( gMC->IsTrackExiting() ||
 474            gMC->IsTrackStop() ||
 475            gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
 476
 477     gMC->TrackPosition(fTrackPosition);
 478     Float_t theta = fTrackMomentum.Theta();
 479     Float_t phi   = fTrackMomentum.Phi();
 480
 481     TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
 482                                fStepSum[idvol]/2.*sin(theta)*sin(phi),
 483                                fStepSum[idvol]/2.*cos(theta),0.0       );
 484     //     if (GetDebug())
 485     //       Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 486     //     if (GetDebug())
 487     //        Info("StepManager Exit ","Track backToWire %f %f %f",backToWire.X(),backToWire.Y(),backToWire.Z()) ;
 488     fTrackPosition-=backToWire;
 489
 490     //-------------- Angle effect
 491     // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
 492
 493     Float_t betaxGamma    = fTrackMomentum.P()/mass;//  pc/mc2
 494     Float_t sigmaEffect10degrees;
 495     Float_t sigmaEffectThetadegrees;
 496     Float_t eLossParticleELossMip;
 497     Float_t yAngleEffect=0.;
 498     Float_t thetawires      =  TMath::Abs( TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) );// We use Pi-theta because z is negative
 499
 500
 501     if (fAngleEffect){
 502     if ( (betaxGamma >3.2)   &&  (thetawires*kRaddeg<=15.) ) {
 503       betaxGamma=TMath::Log(betaxGamma);
 504       eLossParticleELossMip = fElossRatio->Eval(betaxGamma);
 505       // 10 degrees is a reference for a model (arbitrary)
 506       sigmaEffect10degrees=fAngleEffect10->Eval(eLossParticleELossMip);// in micrometers
 507       // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
 508       sigmaEffectThetadegrees =  sigmaEffect10degrees/fAngleEffectNorma->Eval(thetawires*kRaddeg);  // For 5mm gap
 509       if ( (iChamber==1)  ||  (iChamber==2) )
 510         sigmaEffectThetadegrees/=(1.09833e+00+1.70000e-02*(thetawires*kRaddeg)); // The gap is different (4mm)
 511       yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
 512     }
 513     }
 514
 515     // Detection elements ids
 516     AliMUONGeometryModule* geometry
 517       = Chamber(iChamber-1).GetGeometry();
 518
 519     AliMUONGeometryDetElement* detElement
 520       = geometry->FindBySensitiveVolume(CurrentVolumePath());
 521
 522     Int_t detElemId = 0;
 523     if (detElement) detElemId = detElement->GetUniqueID();
 524
 525     if (!detElemId) {
 526       cerr << "Chamber id: "
 527            << setw(3) << iChamber << "  "
 528            << "Current SV: "
 529            <<  CurrentVolumePath()
 530            << "  detElemId: "
 531            << setw(5) << detElemId
 532            << endl;
 533       Double_t x, y, z;
 534       gMC->TrackPosition(x, y, z);
 535       cerr << "   global position: "
 536            << x << ", " << y << ", " << z
 537            << endl;
 538       AliError("DetElemId not identified.");
 539     }
 540
 541     // One hit per chamber
 542     GetMUONData()->AddHit2(fIshunt,
 543                           gAlice->GetMCApp()->GetCurrentTrackNumber(),
 544                           detElemId, ipart,
 545                           fTrackPosition.X(),
 546                           fTrackPosition.Y()+yAngleEffect,
 547                           fTrackPosition.Z(),
 548                           gMC->TrackTime(),
 549                           fTrackMomentum.P(),
 550                           theta,
 551                           phi,
 552                           fStepSum[idvol],
 553                           fDestepSum[idvol],
 554                           fTrackPosition.X(),
 555                           fTrackPosition.Y(),
 556                           fTrackPosition.Z());
 557
 558 //     if (GetDebug()){
 559 //       Info("StepManager Exit","Particle exiting from chamber %d",iChamber);
 560 //       Info("StepManager Exit","StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]);
 561 //       Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 562 //     }
 563     fStepSum[idvol]  =0; // Reset for the next event
 564     fDestepSum[idvol]=0; // Reset for the next event
 565   }
 566 }
 567
 568 //__________________________________________
 569 void AliMUONv1::StepManagerOld()
 570 {
 571   // Old Stepmanager for the chambers
 572   // TBR
 573
 574   Int_t          copy, id;
 575   static Int_t   idvol =-1;
 576   static Int_t   vol[2];
 577   Int_t          ipart;
 578   TLorentzVector pos;
 579   TLorentzVector mom;
 580   Float_t        theta,phi;
 581   Float_t        destep, step;
 582
 583   static Float_t sstep;
 584   static Float_t eloss, eloss2, xhit, yhit, zhit, tof, tlength;
 585   const  Float_t kBig = 1.e10;
 586   static Float_t hits[15];
 587
 588   TClonesArray &lhits = *fHits;
 589
 590   //
 591   //
 592   // Only charged tracks
 593   if( !(gMC->TrackCharge()) ) return;
 594   //
 595   // Only gas gap inside chamber
 596   // Tag chambers and record hits when track enters
 597   id=gMC->CurrentVolID(copy);
 598   vol[0] = GetChamberId(id);
 599   idvol = vol[0] -1;
 600
 601   if (idvol == -1) return;
 602
 603   //
 604   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 605   gMC->TrackPosition(pos);
 606   gMC->TrackMomentum(mom);
 607
 608   ipart  = gMC->TrackPid();
 609
 610   //
 611   // momentum loss and steplength in last step
 612   destep = gMC->Edep();
 613   step   = gMC->TrackStep();
 614   // cout<<"------------"<<step<<endl;
 615   //
 616   // record hits when track enters ...
 617   if( gMC->IsTrackEntering()) {
 618
 619       gMC->SetMaxStep(fMaxStepGas);
 620       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 621       Double_t rt = TMath::Sqrt(tc);
 622       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
 623       Double_t tx = mom[0]/pmom;
 624       Double_t ty = mom[1]/pmom;
 625       Double_t tz = mom[2]/pmom;
 626       Double_t s  = ((AliMUONChamber*)(*fChambers)[idvol])
 627           ->ResponseModel()
 628           ->Pitch()/tz;
 629       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 630       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 631       hits[0] = Float_t(ipart);         // Geant3 particle type
 632       hits[1] = pos[0]+s*tx;            // X-position for hit
 633       hits[2] = pos[1]+s*ty;            // Y-position for hit
 634       hits[3] = pos[2]+s*tz;            // Z-position for hit
 635       hits[4] = theta;                  // theta angle of incidence
 636       hits[5] = phi;                    // phi angle of incidence
 637       hits[8] = 0;//PadHits does not exist anymore  (Float_t) fNPadHits;    // first padhit
 638       hits[9] = -1;                     // last pad hit
 639       hits[10] = mom[3];                // hit momentum P
 640       hits[11] = mom[0];                // Px
 641       hits[12] = mom[1];                // Py
 642       hits[13] = mom[2];                // Pz
 643       tof=gMC->TrackTime();
 644       hits[14] = tof;                   // Time of flight
 645       tlength  = 0;
 646       eloss    = 0;
 647       eloss2   = 0;
 648       sstep=0;
 649       xhit     = pos[0];
 650       yhit     = pos[1];
 651       zhit     = pos[2];
 652       Chamber(idvol).ChargeCorrelationInit();
 653       // Only if not trigger chamber
 654
 655 //       printf("---------------------------\n");
 656 //       printf(">>>> Y =  %f \n",hits[2]);
 657 //       printf("---------------------------\n");
 658
 659
 660
 661      //  if(idvol < AliMUONConstants::NTrackingCh()) {
 662 //        //
 663 //        //  Initialize hit position (cursor) in the segmentation model
 664 //        ((AliMUONChamber*) (*fChambers)[idvol])
 665 //            ->SigGenInit(pos[0], pos[1], pos[2]);
 666 //       } else {
 667 //        //geant3->Gpcxyz();
 668 //        //printf("In the Trigger Chamber #%d\n",idvol-9);
 669 //       }
 670   }
 671   eloss2+=destep;
 672   sstep+=step;
 673
 674   // cout<<sstep<<endl;
 675
 676   //
 677   // Calculate the charge induced on a pad (disintegration) in case
 678   //
 679   // Mip left chamber ...
 680   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 681       gMC->SetMaxStep(kBig);
 682       eloss   += destep;
 683       tlength += step;
 684
 685       Float_t x0,y0,z0;
 686       Float_t localPos[3];
 687       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
 688       gMC->Gmtod(globalPos,localPos,1);
 689
 690       if(idvol < AliMUONConstants::NTrackingCh()) {
 691 // tracking chambers
 692           x0 = 0.5*(xhit+pos[0]);
 693           y0 = 0.5*(yhit+pos[1]);
 694           z0 = 0.5*(zhit+pos[2]);
 695       } else {
 696 // trigger chambers
 697           x0 = xhit;
 698           y0 = yhit;
 699           z0 = 0.;
 700       }
 701
 702
 703       //      if (eloss >0)  MakePadHits(x0,y0,z0,eloss,tof,idvol);
 704
 705
 706       hits[6] = tlength;   // track length
 707       hits[7] = eloss2;    // de/dx energy loss
 708
 709
 710       //      if (fNPadHits > (Int_t)hits[8]) {
 711       //          hits[8] = hits[8]+1;
 712       //          hits[9] = 0: // PadHits does not exist anymore (Float_t) fNPadHits;
 713       //}
 714 //
 715 //    new hit
 716
 717       new(lhits[fNhits++])
 718         AliMUONHit(fIshunt, gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits, kTRUE);
 719       eloss = 0;
 720       //
 721       // Check additional signal generation conditions
 722       // defined by the segmentation
 723       // model (boundary crossing conditions)
 724       // only for tracking chambers
 725   } else if
 726       ((idvol < AliMUONConstants::NTrackingCh()) &&
 727        ((AliMUONChamber*) (*fChambers)[idvol])->SigGenCond(pos[0], pos[1], pos[2]))
 728   {
 729       ((AliMUONChamber*) (*fChambers)[idvol])
 730           ->SigGenInit(pos[0], pos[1], pos[2]);
 731
 732       Float_t localPos[3];
 733       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
 734       gMC->Gmtod(globalPos,localPos,1);
 735
 736       eloss    += destep;
 737
 738       // if (eloss > 0 && idvol < AliMUONConstants::NTrackingCh())
 739       //        MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),pos[2],eloss,tof,idvol);
 740       xhit     = pos[0];
 741       yhit     = pos[1];
 742       zhit     = pos[2];
 743       eloss = 0;
 744       tlength += step ;
 745       //
 746       // nothing special  happened, add up energy loss
 747   } else {
 748       eloss   += destep;
 749       tlength += step ;
 750   }
 751 }
 752
 753 //__________________________________________
 754 void AliMUONv1::StepManagerOld2()
 755 {
 756   // Old Stepmanager for the chambers
 757   Int_t          copy, id;
 758   static Int_t   idvol =-1;
 759   static Int_t   vol[2];
 760   Int_t          ipart;
 761   TLorentzVector pos;
 762   TLorentzVector mom;
 763   Float_t        theta,phi;
 764   Float_t        destep, step;
 765
 766   static Float_t sstep;
 767   static Float_t eloss, eloss2, xhit, yhit, zhit, tof, tlength;
 768   const  Float_t kBig = 1.e10;
 769   static Float_t hits[15];
 770
 771   TClonesArray &lhits = *fHits;
 772
 773   //
 774   //
 775   // Only charged tracks
 776   if( !(gMC->TrackCharge()) ) return;
 777   //
 778   // Only gas gap inside chamber
 779   // Tag chambers and record hits when track enters
 780   id=gMC->CurrentVolID(copy);
 781
 782   Int_t iChamber = GetChamberId(id);
 783   idvol =  iChamber-1;
 784
 785   if (idvol == -1) return;
 786
 787   // Detection elements id
 788   AliMUONGeometryModule* geometry
 789     = Chamber(iChamber-1).GetGeometry();
 790
 791   AliMUONGeometryDetElement* detElement
 792     = geometry->FindBySensitiveVolume(CurrentVolumePath());
 793
 794   Int_t detElemId = 0;
 795   if (detElement) detElemId = detElement->GetUniqueID();
 796
 797   if (!detElemId) {
 798     cerr << "Chamber id: "
 799          << setw(3) << iChamber << "  "
 800          << "Current SV: "
 801          <<  CurrentVolumePath()
 802          << "  detElemId: "
 803          << setw(5) << detElemId
 804          << endl;
 805     Double_t x, y, z;
 806     gMC->TrackPosition(x, y, z);
 807     cerr << "   global position: "
 808          << x << ", " << y << ", " << z
 809          << endl;
 810     AliError("DetElemId not identified.");
 811   }
 812   vol[0] = detElemId;
 813
 814   //
 815   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 816   gMC->TrackPosition(pos);
 817   gMC->TrackMomentum(mom);
 818
 819   ipart  = gMC->TrackPid();
 820
 821   //
 822   // momentum loss and steplength in last step
 823   destep = gMC->Edep();
 824   step   = gMC->TrackStep();
 825   // cout<<"------------"<<step<<endl;
 826   //
 827   // record hits when track enters ...
 828   if( gMC->IsTrackEntering()) {
 829
 830       gMC->SetMaxStep(fMaxStepGas);
 831       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 832       Double_t rt = TMath::Sqrt(tc);
 833       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
 834       Double_t tx = mom[0]/pmom;
 835       Double_t ty = mom[1]/pmom;
 836       Double_t tz = mom[2]/pmom;
 837       Double_t s  = ((AliMUONChamber*)(*fChambers)[idvol])
 838           ->ResponseModel()
 839           ->Pitch()/tz;
 840       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 841       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 842       hits[0] = Float_t(ipart);         // Geant3 particle type
 843       hits[1] = pos[0]+s*tx;            // X-position for hit
 844       hits[2] = pos[1]+s*ty;            // Y-position for hit
 845       hits[3] = pos[2]+s*tz;            // Z-position for hit
 846       hits[4] = theta;                  // theta angle of incidence
 847       hits[5] = phi;                    // phi angle of incidence
 848       hits[8] = 0;//PadHits does not exist anymore  (Float_t) fNPadHits;    // first padhit
 849       hits[9] = -1;                     // last pad hit
 850       hits[10] = mom[3];                // hit momentum P
 851       hits[11] = mom[0];                // Px
 852       hits[12] = mom[1];                // Py
 853       hits[13] = mom[2];                // Pz
 854       tof=gMC->TrackTime();
 855       hits[14] = tof;                   // Time of flight
 856       tlength  = 0;
 857       eloss    = 0;
 858       eloss2   = 0;
 859       sstep=0;
 860       xhit     = pos[0];
 861       yhit     = pos[1];
 862       zhit     = pos[2];
 863       Chamber(idvol).ChargeCorrelationInit();
 864       // Only if not trigger chamber
 865
 866 //       printf("---------------------------\n");
 867 //       printf(">>>> Y =  %f \n",hits[2]);
 868 //       printf("---------------------------\n");
 869
 870
 871
 872      //  if(idvol < AliMUONConstants::NTrackingCh()) {
 873 //        //
 874 //        //  Initialize hit position (cursor) in the segmentation model
 875 //        ((AliMUONChamber*) (*fChambers)[idvol])
 876 //            ->SigGenInit(pos[0], pos[1], pos[2]);
 877 //       } else {
 878 //        //geant3->Gpcxyz();
 879 //        //printf("In the Trigger Chamber #%d\n",idvol-9);
 880 //       }
 881   }
 882   eloss2+=destep;
 883   sstep+=step;
 884
 885   // cout<<sstep<<endl;
 886
 887   //
 888   // Calculate the charge induced on a pad (disintegration) in case
 889   //
 890   // Mip left chamber ...
 891   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 892       gMC->SetMaxStep(kBig);
 893       eloss   += destep;
 894       tlength += step;
 895
 896       Float_t x0,y0,z0;
 897       Float_t localPos[3];
 898       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
 899       gMC->Gmtod(globalPos,localPos,1);
 900
 901       if(idvol < AliMUONConstants::NTrackingCh()) {
 902 // tracking chambers
 903           x0 = 0.5*(xhit+pos[0]);
 904           y0 = 0.5*(yhit+pos[1]);
 905           z0 = 0.5*(zhit+pos[2]);
 906       } else {
 907 // trigger chambers
 908           x0 = xhit;
 909           y0 = yhit;
 910           z0 = 0.;
 911       }
 912
 913
 914       //      if (eloss >0)  MakePadHits(x0,y0,z0,eloss,tof,idvol);
 915
 916
 917       hits[6] = tlength;   // track length
 918       hits[7] = eloss2;    // de/dx energy loss
 919
 920
 921       //      if (fNPadHits > (Int_t)hits[8]) {
 922       //          hits[8] = hits[8]+1;
 923       //          hits[9] = 0: // PadHits does not exist anymore (Float_t) fNPadHits;
 924       //}
 925 //
 926 //    new hit
 927
 928       new(lhits[fNhits++])
 929           AliMUONHit(fIshunt, gAlice->GetMCApp()->GetCurrentTrackNumber(), vol,hits);
 930       eloss = 0;
 931       //
 932       // Check additional signal generation conditions
 933       // defined by the segmentation
 934       // model (boundary crossing conditions)
 935       // only for tracking chambers
 936   } else if
 937       ((idvol < AliMUONConstants::NTrackingCh()) &&
 938        ((AliMUONChamber*) (*fChambers)[idvol])->SigGenCond(pos[0], pos[1], pos[2]))
 939   {
 940       ((AliMUONChamber*) (*fChambers)[idvol])
 941           ->SigGenInit(pos[0], pos[1], pos[2]);
 942
 943       Float_t localPos[3];
 944       Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
 945       gMC->Gmtod(globalPos,localPos,1);
 946
 947       eloss    += destep;
 948
 949       // if (eloss > 0 && idvol < AliMUONConstants::NTrackingCh())
 950       //        MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),pos[2],eloss,tof,idvol);
 951       xhit     = pos[0];
 952       yhit     = pos[1];
 953       zhit     = pos[2];
 954       eloss = 0;
 955       tlength += step ;
 956       //
 957       // nothing special  happened, add up energy loss
 958   } else {
 959       eloss   += destep;
 960       tlength += step ;
 961   }
 962 }