#include <TRandom.h>
#include <TTUBE.h>
#include <TVirtualMC.h>
+#include <TParticle.h>
#include "AliCallf77.h"
#include "AliConst.h"
//___________________________________________
AliMUONv1::AliMUONv1() : AliMUON()
+ ,fTrackMomentum(), fTrackPosition()
{
// Constructor
- fChambers = 0;
- fStations = 0;
- fStepManagerVersionOld = kFALSE;
- fStepManagerVersionNew = kFALSE;
- fStepManagerVersionTest = kFALSE;
-
- fStepMaxInActiveGas = 2.0;
-}
-
-
+ fChambers = 0;
+ fStations = 0;
+ fStepManagerVersionOld = kFALSE;
+ fStepMaxInActiveGas = 0.6;
+ fStepSum = 0x0;
+ fDestepSum = 0x0;
+ fElossRatio = 0x0;
+ fAngleEffect10 = 0x0;
+ fAngleEffectNorma= 0x0;
+}
//___________________________________________
AliMUONv1::AliMUONv1(const char *name, const char *title)
- : AliMUON(name,title)
+ : AliMUON(name,title), fTrackMomentum(), fTrackPosition()
{
// Constructor
// By default include all stations
factory.Build(this, title);
fStepManagerVersionOld = kFALSE;
- fStepManagerVersionNew = kFALSE;
- fStepManagerVersionTest = kFALSE;
- fStepMaxInActiveGas = 2.0;
+ fStepMaxInActiveGas = 0.6;
+
+ fStepSum = new Float_t [AliMUONConstants::NCh()];
+ fDestepSum = new Float_t [AliMUONConstants::NCh()];
+ for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
+ fStepSum[i] =0.0;
+ fDestepSum[i]=0.0;
+ }
+ // Ratio of particle mean eloss with respect MIP's Khalil Boudjemline, sep 2003, PhD.Thesis and Particle Data Book
+ fElossRatio = new TF1("ElossRatio","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",0.5,5.);
+ fElossRatio->SetParameter(0,1.02138);
+ fElossRatio->SetParameter(1,-9.54149e-02);
+ fElossRatio->SetParameter(2,+7.83433e-02);
+ fElossRatio->SetParameter(3,-9.98208e-03);
+ fElossRatio->SetParameter(4,+3.83279e-04);
+
+ // Angle effect in tracking chambers at theta =10 degres as a function of ElossRatio (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) (in micrometers)
+ fAngleEffect10 = new TF1("AngleEffect10","[0]+[1]*x+[2]*x*x",0.5,3.0);
+ fAngleEffect10->SetParameter(0, 1.90691e+02);
+ fAngleEffect10->SetParameter(1,-6.62258e+01);
+ fAngleEffect10->SetParameter(2,+1.28247e+01);
+ // Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis)
+ // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
+ fAngleEffectNorma = new TF1("AngleEffectNorma","[0]+[1]*x+[2]*x*x+[3]*x*x*x",0.0,10.0);
+ fAngleEffectNorma->SetParameter(0,4.148);
+ fAngleEffectNorma->SetParameter(1,-6.809e-01);
+ fAngleEffectNorma->SetParameter(2,5.151e-02);
+ fAngleEffectNorma->SetParameter(3,-1.490e-03);
}
//___________________________________________
//cp
}
-//___________________________________________
+
+//_______________________________________________________________________________
+Int_t AliMUONv1::GetChamberId(Int_t volId) const
+{
+// Check if the volume with specified volId is a sensitive volume (gas)
+// of some chamber and returns the chamber number;
+// if not sensitive volume - return 0.
+// ---
+
+ for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
+ if (volId==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) return i;
+
+ return 0;
+}
+//_______________________________________________________________________________
void AliMUONv1::StepManager()
{
- if (fStepManagerVersionOld) {
+ if (fStepManagerVersionOld) {
StepManagerOld();
return;
}
- if (fStepManagerVersionNew) {
- StepManagerNew();
- return;
- }
-
- if (fStepManagerVersionTest) {
- StepManagerTest();
- return;
- }
-
-
- // Volume id
- Int_t copy, id;
- Int_t idvol;
- Int_t iChamber=0;
- // Particule id, pos and mom vectors,
- // theta, phi angles with respect the normal of the chamber,
- // spatial step, delta_energy and time of flight
- Int_t ipart;
- TLorentzVector pos, mom;
- Float_t theta, phi, tof;
- Float_t destep, step;
- const Float_t kBig = 1.e10;
// Only charged tracks
if( !(gMC->TrackCharge()) ) return;
-
+ // Only charged tracks
+
// Only gas gap inside chamber
// Tag chambers and record hits when track enters
- idvol=-1;
+ Int_t idvol=-1;
+ Int_t iChamber=0;
+ Int_t id=0;
+ Int_t copy;
+ const Float_t kBig = 1.e10;
+
id=gMC->CurrentVolID(copy);
+ // printf("id == %d \n",id);
for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++) {
if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) {
iChamber = i;
idvol = i-1;
}
}
- if (idvol == -1) return;
-
- // printf(">>>> This Chamber %d\n",iChamber);
-
- // record hits when track enters ...
- if( gMC->IsTrackEntering()) gMC->SetMaxStep(fStepMaxInActiveGas);
-
- if (gMC->TrackStep() > 0.) {
- // Get current particle id (ipart), track position (pos) and momentum (mom)
- gMC->TrackPosition(pos);
- gMC->TrackMomentum(mom);
- ipart = gMC->TrackPid();
- theta = mom.Theta()*kRaddeg; // theta of track
- phi = mom.Phi() *kRaddeg; // phi of the track
- tof = gMC->TrackTime(); // Time of flight
- //
- // momentum loss and steplength in last step
- destep = gMC->Edep();
- step = gMC->TrackStep();
-
- //new hit
- GetMUONData()->AddHit(fIshunt, gAlice->GetCurrentTrackNumber(), iChamber, ipart,
- pos.X(), pos.Y(), pos.Z(), tof, mom.P(),
- theta, phi, step, destep);
- }
- // Track left chamber ...
- if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
- gMC->SetMaxStep(kBig);
+ if (idvol == -1) {
+ return;
}
-}
-
-
-Int_t AliMUONv1::GetChamberId(Int_t volId) const
-{
-// Check if the volume with specified volId is a sensitive volume (gas)
-// of some chamber and returns the chamber number;
-// if not sensitive volume - return 0.
-// ---
- for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
- if (volId==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) return i;
-
- return 0;
-}
-//__
-
-
-
-void AliMUONv1::StepManagerTest()
-{
- return;
-}
-//________________________________________
-void AliMUONv1::StepManagerNew()
-{
-
-// // Volume id
-// Int_t copy, id;
-// Int_t idvol;
-// Int_t iChamber=0;
-// // Particule id, pos and mom vectors,
-// // theta, phi angles with respect the normal of the chamber,
-// // spatial step, delta_energy and time of flight
-// Int_t ipart;
-// TLorentzVector pos, mom;
-// Float_t theta, phi, tof;
-// Float_t destep, step;
-// const Float_t kBig = 1.e10;
-
-// // Only charged tracks
-// if( !(gMC->TrackCharge()) ) return;
-
-// // Only gas gap inside chamber
-// // Tag chambers and record hits when track enters
-// idvol=-1;
-// id=gMC->CurrentVolID(copy);
-// for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++) {
-// if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) {
-// iChamber = i;
-// idvol = i-1;
-// }
-// }
-// static Float_t Sstep[20]; // Sum of steps per chamber
-// // static Float_t Sdestep[20]; // Sum of eloss per chamber
-// Float_t GAP;
-// Float_t TEST;
-
-// if (idvol == -1) return;
-
-// // printf(">>>> This Chamber %d\n",iChamber);
-
-// // record hits when track enters ...
-// //if( gMC->IsTrackEntering()) gMC->SetMaxStep(fStepMaxInActiveGas);
-
-// if (gMC->TrackStep() > 0.) {
-// // Get current particle id (ipart), track position (pos) and momentum (mom)
-// gMC->TrackPosition(pos);
-// gMC->TrackMomentum(mom);
-// ipart = gMC->TrackPid(); // Particle
-// theta = mom.Theta()*kRaddeg; // theta of track
-// phi = mom.Phi() *kRaddeg; // phi of the track
-// tof = gMC->TrackTime(); // Time of flight
-// //
-// // momentum loss and steplength in last step
-// destep = gMC->Edep();
-// step = gMC->TrackStep();
-
-// Sstep[iChamber]+=step;
-// // Sdestep[iChamber]+=destep;
+ if( gMC->IsTrackEntering() ) {
+ Float_t theta = fTrackMomentum.Theta();
+ Float_t phi = fTrackMomentum.Phi();
+ Float_t theta_wires = TMath::Abs( TMath::ASin( TMath::Sin(theta) * TMath::Sin(phi) ) );
+ if ( (theta_wires*kRaddeg>=10) ) gMC->SetMaxStep(fStepMaxInActiveGas);
+ }
-// }
+// if (GetDebug()) {
+// Float_t z = ( (AliMUONChamber*)(*fChambers)[idvol])->Z() ;
+// Info("StepManager Step","Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber, z);
+// }
+ // Particule id and mass,
+ Int_t ipart = gMC->TrackPid();
+ Float_t mass = gMC->TrackMass();
+
+ fDestepSum[idvol]+=gMC->Edep();
+ // Get current particle id (ipart), track position (pos) and momentum (mom)
+ if ( fStepSum[idvol]==0.0 ) gMC->TrackMomentum(fTrackMomentum);
+ fStepSum[idvol]+=gMC->TrackStep();
-// step = Sstep[iChamber]; // Total step >= gap
-// // destep = Sdestep[iChamber]; // Total eloss
-
-
-// // Track left chamber ...
-// if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
-// gMC->SetMaxStep(kBig);
-
-// Sstep[iChamber]=0; // Reset for the next event
-// //Sdestep[iChamber]=0; // Reset for the next event
+// if (GetDebug()) {
+// Info("StepManager Step","iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g",
+// 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()) ;
+// gMC->TrackPosition(fTrackPosition);
+// Info("StepManager Step","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
+// }
-// if (iChamber>=1 && iChamber<=2) GAP=0.4;
-// if (iChamber>=11 && iChamber<=14) GAP=0.2;
-// if (iChamber>=3 && iChamber<=10) GAP=0.5;
-
-// TF1 *ELOSS1 = new TF1("Gauss1","exp(-((x-4.13727e+01)**2)/(2*1.42223e+01**2))",0,75);
-// TF1 *ELOSS2 = new TF1("Gauss2","exp(-((x+6.83795e+02)**2)/(2*4.48415e+02**2))",75,350);
-// TEST=gRandom->Rndm();
-// if (TEST <=0.89) destep=ELOSS1->GetRandom();
-// else destep=ELOSS2->GetRandom();
-// destep*=pow(10,-6)*0.0274;
-// destep*=GAP/0.5;
-
-// // One hit per chamber
-// GetMUONData()->AddHit(fIshunt, gAlice->GetCurrentTrackNumber(), iChamber, ipart,
-// pos.X()-(step/2*sin(theta*kDegrad)*cos(phi*kDegrad)), pos.Y()-(step/2*sin(theta*kDegrad)*sin(phi*kDegrad)), pos.Z()-GAP/2, tof, mom.P(),theta, phi, step, destep);
+ // Track left chamber or StepSum larger than fStepMaxInActiveGas
+ if ( gMC->IsTrackExiting() ||
+ gMC->IsTrackStop() ||
+ gMC->IsTrackDisappeared()||
+ (fStepSum[idvol]>fStepMaxInActiveGas) ) {
+
+ if ( gMC->IsTrackExiting() ||
+ gMC->IsTrackStop() ||
+ gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
-// }
+ gMC->TrackPosition(fTrackPosition);
+ Float_t theta = fTrackMomentum.Theta();
+ Float_t phi = fTrackMomentum.Phi();
+
+ TLorentzVector BackToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
+ fStepSum[idvol]/2.*sin(theta)*sin(phi),
+ fStepSum[idvol]/2.*cos(theta),0.0 );
+ // if (GetDebug())
+ // Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
+ // if (GetDebug())
+ // Info("StepManager Exit ","Track BackToWire %f %f %f",BackToWire.X(),BackToWire.Y(),BackToWire.Z()) ;
+ fTrackPosition-=BackToWire;
+
+ //-------------- Angle effect
+ // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
+
+ Float_t Beta_x_Gamma = fTrackMomentum.P()/mass;// pc/mc2
+ Float_t SigmaEffect_10degrees;
+ Float_t SigmaEffect_thetadegrees;
+ Float_t ELossParticle_ELossMip;
+ Float_t YAngleEffect=0.;
+ Float_t theta_wires = TMath::Abs( TMath::ASin( TMath::Sin(theta) * TMath::Sin(phi) ) );
+
+ if ( (Beta_x_Gamma >3.2) && (theta_wires*kRaddeg<=10) ) {
+ Beta_x_Gamma=TMath::Log(Beta_x_Gamma);
+ ELossParticle_ELossMip = fElossRatio->Eval(Beta_x_Gamma);
+ // 10 degrees is a reference for a model (arbitrary)
+ SigmaEffect_10degrees=fAngleEffect10->Eval(ELossParticle_ELossMip);// in micrometers
+ // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
+ SigmaEffect_thetadegrees = SigmaEffect_10degrees/fAngleEffectNorma->Eval(theta_wires*kRaddeg); // For 5mm gap
+ if ( (iChamber==1) || (iChamber==2) )
+ SigmaEffect_thetadegrees/=(1.09833e+00+1.70000e-02*theta_wires*kRaddeg); // The gap is different (4mm)
+ YAngleEffect=1.e-04*gRandom->Gaus(0,SigmaEffect_thetadegrees); // Error due to the angle effect in cm
+ }
+
+
+ // One hit per chamber
+ GetMUONData()->AddHit(fIshunt, gAlice->GetCurrentTrackNumber(), iChamber, ipart,
+ fTrackPosition.X(), fTrackPosition.Y()+YAngleEffect, fTrackPosition.Z(), 0.0,
+ fTrackMomentum.P(),theta, phi, fStepSum[idvol], fDestepSum[idvol],
+ fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z());
+// if (GetDebug()){
+// Info("StepManager Exit","Particle exiting from chamber %d",iChamber);
+// Info("StepManager Exit","StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]);
+// Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
+// }
+ fStepSum[idvol] =0; // Reset for the next event
+ fDestepSum[idvol]=0; // Reset for the next event
+ }
}
//___________________________________________
{
//
// Config file for MUON test
- // Gines MARITNEZ, Subatech, mai 2003, august 2003
+ // Gines MARTINEZ, Subatech, mai 2003, august 2003
//
//=====================================================================
if (!strcmp(option,"box")) {
AliGenBox * gener = new AliGenBox(1);
gener->SetMomentumRange(20.,20.1);
- gener->SetPhiRange(45., 45.01);
- gener->SetThetaRange(4.000,4.001);
+ gener->SetPhiRange(-180, 180.01);
+ gener->SetThetaRange(171.000,178.001);
gener->SetPart(13); // Muons
gener->SetOrigin(0.,0., 0.); //vertex position
gener->SetSigma(0.0, 0.0, 0.0); //Sigma in (X,Y,Z) (cm) on IP position
gener->SetMomentumRange(0,999);
gener->SetPtRange(0,100.);
gener->SetPhiRange(-180, 180);
- gener->SetYRange(2.5,4);
gener->SetCutOnChild(1);
- gener->SetChildThetaRange(2.0,9);
+ gener->SetChildThetaRange(171.0,178.0);
gener->SetOrigin(0,0,0); //vertex position gener->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position
gener->SetForceDecay(kDiMuon);
gener->SetTrackingFlag(1);
//=============================================================
// Field (L3 0.4 T)
- AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k4kG);
- gAlice->SetField(field);
- //gAlice->SetField(2,1) ; //(-999,2);
-
+ AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k4kG);
+ gAlice->SetField(field);
//=================== Alice BODY parameters =============================
AliBODY *BODY = new AliBODY("BODY","Alice envelop");
//=================== MUON Subsystem ===========================
- AliMUONv1 *MUON = new AliMUONv1("MUON","default");
-
+ cout << ">>> Config_MUON_test.C: Creating AliMUONv1 ..."<<endl;
+ AliMUONv1 *MUON = new AliMUONv1("MUON","default");
}
-
Float_t EtaToTheta(Float_t arg){
return (180./TMath::Pi())*2.*atan(exp(-arg));
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff