#include <TLorentzVector.h>
#include <TMath.h>
#include <TRandom.h>
#include <TVirtualMC.h>
#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include <TGeoPhysicalNode.h>
#include "AliTrackReference.h"
#include "AliMC.h"
#include "AliRun.h"
#include "AliGeomManager.h"
#include "AliTRDgeometry.h"
#include "AliTRDCommonParam.h"
#include "AliTRDsimTR.h"
#include "AliTRDtestG4.h"
ClassImp(AliTRDtestG4)
AliTRDtestG4::AliTRDtestG4()
:AliTRD()
,fTRon(kTRUE)
,fTR(NULL)
,fStepSize(0)
,fWion(0)
{
}
AliTRDtestG4::AliTRDtestG4(const char *name, const char *title)
:AliTRD(name,title)
,fTRon(kTRUE)
,fTR(NULL)
,fStepSize(0.1)
,fWion(0)
{
SetBufferSize(128000);
if (AliTRDCommonParam::Instance()->IsXenon()) {
fWion = 23.53;
}
else if (AliTRDCommonParam::Instance()->IsArgon()) {
fWion = 27.21;
}
else {
AliFatal("Wrong gas mixture");
exit(1);
}
}
AliTRDtestG4::~AliTRDtestG4()
{
if (fTR) {
delete fTR;
fTR = 0;
}
}
void AliTRDtestG4::AddAlignableVolumes() const
{
TString volPath;
TString symName;
TString vpStr = "ALIC_1/B077_1/BSEGMO";
TString vpApp1 = "_1/BTRD";
TString vpApp2 = "_1";
TString vpApp3a = "/UTR1_1/UTS1_1/UTI1_1/UT";
TString vpApp3b = "/UTR2_1/UTS2_1/UTI2_1/UT";
TString vpApp3c = "/UTR3_1/UTS3_1/UTI3_1/UT";
TString vpApp3d = "/UTR4_1/UTS4_1/UTI4_1/UT";
TString snStr = "TRD/sm";
TString snApp1 = "/st";
TString snApp2 = "/pl";
for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
volPath = vpStr;
volPath += isector;
volPath += vpApp1;
volPath += isector;
volPath += vpApp2;
symName = snStr;
symName += Form("%02d",isector);
gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
}
AliGeomManager::ELayerID idTRD1 = AliGeomManager::kTRD1;
Int_t layer, modUID;
for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
if (fGeometry->GetSMstatus(isector) == 0) continue;
for (Int_t istack = 0; istack < AliTRDgeometry::Nstack(); istack++) {
for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) {
layer = idTRD1 + ilayer;
modUID = AliGeomManager::LayerToVolUIDSafe(layer,isector*5+istack);
Int_t idet = AliTRDgeometry::GetDetectorSec(ilayer,istack);
volPath = vpStr;
volPath += isector;
volPath += vpApp1;
volPath += isector;
volPath += vpApp2;
switch (isector) {
case 17:
if ((istack == 4) && (ilayer == 4)) {
continue;
}
volPath += vpApp3d;
break;
case 13:
case 14:
case 15:
if (istack == 2) {
continue;
}
volPath += vpApp3c;
break;
case 11:
case 12:
volPath += vpApp3b;
break;
default:
volPath += vpApp3a;
};
volPath += Form("%02d",idet);
volPath += vpApp2;
symName = snStr;
symName += Form("%02d",isector);
symName += snApp1;
symName += istack;
symName += snApp2;
symName += ilayer;
TGeoPNEntry *alignableEntry =
gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data(),modUID);
if (alignableEntry) {
TGeoHMatrix *globMatrix = alignableEntry->GetGlobalOrig();
Double_t sectorAngle = 20.0 * (isector % 18) + 10.0;
TGeoHMatrix *t2lMatrix = new TGeoHMatrix();
t2lMatrix->RotateZ(sectorAngle);
t2lMatrix->MultiplyLeft(&(globMatrix->Inverse()));
alignableEntry->SetMatrix(t2lMatrix);
}
else {
AliError(Form("Alignable entry %s is not valid!",symName.Data()));
}
}
}
}
}
void AliTRDtestG4::CreateGeometry()
{
AliModule* frame = gAlice->GetModule("FRAME");
if (!frame) {
AliError("TRD needs FRAME to be present\n");
return;
}
AliTRD::CreateGeometry();
}
void AliTRDtestG4::CreateMaterials()
{
AliTRD::CreateMaterials();
}
void AliTRDtestG4::CreateTRhit(Int_t det)
{
const Int_t kNTR = 50;
TLorentzVector mom;
TLorentzVector pos;
Float_t eTR[kNTR];
Int_t nTR;
TVirtualMC::GetMC()->TrackMomentum(mom);
Float_t pTot = mom.Rho();
fTR->CreatePhotons(11,pTot,nTR,eTR);
if (nTR > kNTR) {
AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
}
for (Int_t iTR = 0; iTR < nTR; iTR++) {
Float_t energyMeV = eTR[iTR] * 0.001;
Float_t energyeV = eTR[iTR] * 1000.0;
Float_t absLength = 0.0;
Float_t sigma = 0.0;
Double_t muMy = fTR->GetMuMy(energyMeV);
sigma = muMy * fFoilDensity;
if (sigma > 0.0) {
absLength = gRandom->Exp(1.0/sigma);
if (absLength < AliTRDgeometry::MyThick()) {
continue;
}
}
else {
continue;
}
Double_t muNo = 0.0;
if (AliTRDCommonParam::Instance()->IsXenon()) {
muNo = fTR->GetMuXe(energyMeV);
}
else if (AliTRDCommonParam::Instance()->IsArgon()) {
muNo = fTR->GetMuAr(energyMeV);
}
Double_t muCO = fTR->GetMuCO(energyMeV);
sigma = (fGasNobleFraction * muNo + (1.0 - fGasNobleFraction) * muCO)
* fGasDensity
* fTR->GetTemp();
if (sigma > 0.0) {
absLength = gRandom->Exp(1.0/sigma);
if (absLength > (AliTRDgeometry::DrThick()
+ AliTRDgeometry::AmThick())) {
continue;
}
}
else {
continue;
}
Float_t posHit[3];
TVirtualMC::GetMC()->TrackPosition(pos);
posHit[0] = pos[0] + mom[0] / pTot * absLength;
posHit[1] = pos[1] + mom[1] / pTot * absLength;
posHit[2] = pos[2] + mom[2] / pTot * absLength;
Int_t q = ((Int_t) (energyeV / fWion));
AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
,det
,posHit
,-q
,TVirtualMC::GetMC()->TrackTime()*1.0e06
,kTRUE);
}
}
void AliTRDtestG4::Init()
{
AliTRD::Init();
AliDebug(1,"Slow simulator\n");
if (!fTRon) {
AliInfo("TR simulation off");
}
else {
fTR = new AliTRDsimTR();
}
AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
}
void AliTRDtestG4::StepManager()
{
const Int_t kPdgElectron = 11;
Int_t layer = 0;
Int_t stack = 0;
Int_t sector = 0;
Int_t det = 0;
Int_t qTot;
Float_t hits[3];
Double_t eDep;
Bool_t drRegion = kFALSE;
Bool_t amRegion = kFALSE;
TString cIdPath;
Char_t cIdSector[3];
cIdSector[2] = 0;
TString cIdCurrent;
TString cIdSensDr = "J";
TString cIdSensAm = "K";
Char_t cIdChamber[3];
cIdChamber[2] = 0;
TLorentzVector pos;
TLorentzVector mom;
const Int_t kNlayer = AliTRDgeometry::Nlayer();
const Int_t kNstack = AliTRDgeometry::Nstack();
const Int_t kNdetsec = kNlayer * kNstack;
const Double_t kBig = 1.0e+12;
const Float_t kEkinMinStep = 1.0e-5;
const Double_t kScaleG4 = 1.12;
if (!fPrimaryIonisation) TVirtualMC::GetMC()->SetMaxStep(kBig);
if ((!TVirtualMC::GetMC()->TrackCharge()) ||
TVirtualMC::GetMC()->IsTrackDisappeared()) {
return;
}
cIdCurrent = TVirtualMC::GetMC()->CurrentVolName();
if (cIdSensDr == cIdCurrent[1]) {
drRegion = kTRUE;
}
if (cIdSensAm == cIdCurrent[1]) {
amRegion = kTRUE;
}
if ((!drRegion) &&
(!amRegion)) {
return;
}
TVirtualMC::GetMC()->TrackPosition(pos);
hits[0] = pos[0];
hits[1] = pos[1];
hits[2] = pos[2];
cIdPath = gGeoManager->GetPath();
cIdSector[0] = cIdPath[21];
cIdSector[1] = cIdPath[22];
sector = atoi(cIdSector);
cIdChamber[0] = cIdCurrent[2];
cIdChamber[1] = cIdCurrent[3];
Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
stack = ((Int_t) idChamber / kNlayer);
layer = ((Int_t) idChamber % kNlayer);
det = fGeometry->GetDetector(layer,stack,sector);
Int_t trkStat = 0;
if ((drRegion) &&
(TVirtualMC::GetMC()->IsTrackEntering())) {
TVirtualMC::GetMC()->TrackMomentum(mom);
AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
trkStat = 1;
if ((fTR) &&
(fTRon) &&
(TMath::Abs(TVirtualMC::GetMC()->TrackPid()) == kPdgElectron)) {
CreateTRhit(det);
}
}
else if ((amRegion) &&
(TVirtualMC::GetMC()->IsTrackExiting())) {
TVirtualMC::GetMC()->TrackMomentum(mom);
AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
trkStat = 2;
}
eDep = TMath::Max(TVirtualMC::GetMC()->Edep(),0.0) * 1.0e+09;
eDep /= kScaleG4;
qTot = (Int_t) (eDep / fWion);
if ((qTot) ||
(trkStat)) {
AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
,det
,hits
,qTot
,TVirtualMC::GetMC()->TrackTime()*1.0e06
,drRegion);
}
if ((TVirtualMC::GetMC()->Etot() - TVirtualMC::GetMC()->TrackMass()) < kEkinMinStep) {
return;
}
if (!fPrimaryIonisation) TVirtualMC::GetMC()->SetMaxStep(fStepSize);
}