using std::ifstream;
using std::ios_base;
+extern "C"{
+ Gas gaspar_;
+};
ClassImp(AliTPCv2)
//_____________________________________________________________________________
SetBufferSize(128000);
+ if(!fTPCParam) {AliFatal("TPC parameters not set");
+ return;
+ }
+ gaspar_.fpot=fTPCParam->GetFpot();
+ gaspar_.eend=fTPCParam->GetEend();
+ gaspar_.eexpo=fTPCParam->GetExp();
-// if (fTPCParam)
-// fTPCParam->Write(fTPCParam->GetTitle());
}
//_____________________________________________________________________________
TGeoCombiTrans *transf[13];
Char_t name[30];
for(Int_t i=0;i<13;i++){
- //sprintf(name,"transf%d",i);
snprintf(name,30,"transf%d",i);
transf[i]= new TGeoCombiTrans(name,0.,-2.,-9.+i*1.5,rhole);
transf[i]->RegisterYourself();
// union expression for holes
TString operl("hhole:transf0");
for (Int_t i=1;i<13;i++){
- //sprintf(name,"+hhole:transf%d",i);
snprintf(name,30,"+hhole:transf%d",i);
operl.Append(name);
}
//
TString opers("hhole:transf1");
for (Int_t i=2;i<12;i++){
- //sprintf(name,"+hhole:transf%d",i);
snprintf(name,30,"+hhole:transf%d",i);
opers.Append(name);
}
TGeoRotation *rotpos[2];
//
TGeoRotation *rotrod1[2];
+ //
+ // clamps holding rods
+ //
+ TGeoBBox *clampi1 = new TGeoBBox("clampi1",0.2,3.1,0.8);
+ TGeoVolume *clampi1v = new TGeoVolume("TPC_clampi1v",clampi1,m6);
+ //
+ pointstrap[0]=0.49;
+ pointstrap[1]=0.375;
+ //
+ pointstrap[2]=0.49;
+ pointstrap[3]=-0.375;
+ //
+ pointstrap[4]=-0.49;
+ pointstrap[5]=-0.375;
+ //
+ pointstrap[6]=-0.49;
+ pointstrap[7]=1.225;
+ //
+ pointstrap[8]=0.49;
+ pointstrap[9]=0.375;
+ //
+ pointstrap[10]=0.49;
+ pointstrap[11]=-0.375;
+ //
+ pointstrap[12]=-0.49;
+ pointstrap[13]=-0.375;
+ //
+ pointstrap[14]=-0.49;
+ pointstrap[15]=1.225;
+ //
+ TGeoArb8 *clitrap = new TGeoArb8("clitrap",0.25,pointstrap);
+ TGeoVolume *clitrapv = new TGeoVolume("TPC_clitrapv",clitrap,m6);
+ //
+ TGeoRotation *clamprot = new TGeoRotation();
+ clamprot->RotateX(180.);
+ //
+ new TGeoBBox("clibox",1.125,3.1,.1);
+ new TGeoTube("clitub",0.,2.2,0.1);
+ //
+ // copmisite shape for the clamp holder
+ //
+ TGeoTranslation *clitr1 = new TGeoTranslation("clitr1",1.125,0.,0.);
+ clitr1->RegisterYourself();
+ TGeoCompositeShape *clihold = new TGeoCompositeShape("clihold","clibox-clitub:clitr1");
+ TGeoVolume *cliholdv = new TGeoVolume("TPC_cliholdv",clihold,m6);
+ //
+ // now assembly the whole inner clamp
+ //
+ TGeoVolume *iclamp = new TGeoVolumeAssembly("TPC_iclamp");
+ //
+ iclamp->AddNode(clampi1v,1); //main box
+ iclamp->AddNode(clitrapv,1,new TGeoTranslation(0.69,-2.725,0.35)); //trapezoids
+ iclamp->AddNode(clitrapv,2,new TGeoTranslation(0.69,-2.725,-0.35));
+ iclamp->AddNode(clitrapv,3,new TGeoCombiTrans(0.69,2.725,0.35,clamprot));
+ iclamp->AddNode(clitrapv,4,new TGeoCombiTrans(0.69,2.725,-0.35,clamprot));
+ iclamp->AddNode(cliholdv,1,new TGeoTranslation(1.325,0.,0.)); //holder
+ //
+ // outer clamps
+ //
+ TGeoBBox *clampo1 = new TGeoBBox("clampo1",0.25,3.1,1.);
+ TGeoBBox *clampo2 = new TGeoBBox("clampo2",0.4,0.85,1.);
+ //
+ TGeoVolume *clampo1v = new TGeoVolume("TPC_clampo1v",clampo1,m6);
+ TGeoVolume *clampo2v = new TGeoVolume("TPC_clampo2v",clampo2,m6);
+ //
+ TGeoVolumeAssembly *oclamp = new TGeoVolumeAssembly("TPC_oclamp");
+ //
+ oclamp->AddNode(clampo1v,1);
+ //
+ oclamp->AddNode(clampo2v,1,new TGeoTranslation(0.65,-2.25,0));
+ oclamp->AddNode(clampo2v,2,new TGeoTranslation(0.65,2.25,0));
+
+ //
+ pointstrap[0]=0.375;
+ pointstrap[1]=0.75;
+ pointstrap[2]=0.375;
+ pointstrap[3]=-0.35;
+ pointstrap[5]=-0.375;
+ pointstrap[4]=-0.35;
+ pointstrap[6]=-0.375;
+ pointstrap[7]=0.35;
+ //
+ pointstrap[8]=0.375;
+ pointstrap[9]=0.75;
+ pointstrap[10]=0.375;
+ pointstrap[11]=-0.35;
+ pointstrap[12]=-0.375;
+ pointstrap[13]=-0.35;
+ pointstrap[14]=-0.375;
+ pointstrap[15]=0.35;
+ //
+ TGeoArb8 *clotrap = new TGeoArb8("clotrap",0.25,pointstrap);
+ TGeoVolume *clotrapv = new TGeoVolume("TPC_clotrapv",clotrap,m6);
+ //
+ oclamp->AddNode(clotrapv,1,new TGeoTranslation(-0.625,-2.75,0.35));
+ oclamp->AddNode(clotrapv,2,new TGeoTranslation(-0.625,-2.75,-0.35));
+ oclamp->AddNode(clotrapv,3,new TGeoCombiTrans(-0.625,2.75,0.35,clamprot));
+ oclamp->AddNode(clotrapv,4,new TGeoCombiTrans(-0.625,2.75,-0.35,clamprot));
+ //
+ TGeoBBox *clampo3 = new TGeoBBox("clampo3",1.6,0.45,.1);
+ TGeoVolume *clampo3v = new TGeoVolume("TPC_clampo3v",clampo3,m6);
+ //
+ oclamp->AddNode(clampo3v,1,new TGeoTranslation(-1.85,2.625,0.));
+ oclamp->AddNode(clampo3v,2,new TGeoTranslation(-1.85,-2.625,0));
+ //
+ TGeoTubeSeg *clampo4 = new TGeoTubeSeg("clampo4",2.2,3.1,0.1,90.,270.);
+ TGeoVolume *clampo4v = new TGeoVolume("TPC_clampo4v",clampo4,m6);
+ //
+ oclamp->AddNode(clampo4v,1,new TGeoTranslation(-3.45,0.,0.));
+
+
-
//v9 - drift gas
TGeoRotation rot102("rot102");
v9->AddNode(tpcihpl,i+1,new TGeoCombiTrans(x, y, 0., rot12));
//
if(i==11){//resistor rod inner
- rotrod.RotateZ(-90.+angle);
+ rotrod.RotateZ(-90.+i*20.);
rotrod1[0]= new TGeoRotation();
rotpos[0]= new TGeoRotation();
//
- rotrod1[0]->RotateZ(-90.+angle);
+ rotrod1[0]->RotateZ(90.+i*20.);
*rotpos[0] = refl*rotrod; //rotation+reflection
v9->AddNode(tpcrrod,1,new TGeoCombiTrans(x,y, z, rotrod1[0])); //A
v9->AddNode(tpcrrod,2,new TGeoCombiTrans(x,y,-z, rotpos[0])); //C
v9->AddNode(tpcmrod,i+1,new TGeoTranslation(x,y,z));//shaft
v9->AddNode(tpcmrod,i+19,new TGeoCombiTrans(x,y,-z,ref));//muon
}
+ //
+ // inner clamps positioning
+ //
+ r=79.05;
+ x=r * TMath::Cos(angle);
+ y=r * TMath::Sin(angle);
+ rot12= new TGeoRotation();
+ rot12->RotateZ(i*20.);
+ //
+ //A-side
+ v9->AddNode(iclamp,7*i+1,new TGeoCombiTrans(x,y,5.25,rot12));
+ v9->AddNode(iclamp,7*i+2,new TGeoCombiTrans(x,y,38.25,rot12));
+ v9->AddNode(iclamp,7*i+3,new TGeoCombiTrans(x,y,80.25,rot12));
+ v9->AddNode(iclamp,7*i+4,new TGeoCombiTrans(x,y,122.25,rot12));
+ v9->AddNode(iclamp,7*i+5,new TGeoCombiTrans(x,y,164.25,rot12));
+ v9->AddNode(iclamp,7*i+6,new TGeoCombiTrans(x,y,206.25,rot12));
+ v9->AddNode(iclamp,7*i+7,new TGeoCombiTrans(x,y,246.75,rot12));
+ //C-side
+ v9->AddNode(iclamp,7*i+127,new TGeoCombiTrans(x,y,-5.25,rot12));
+ v9->AddNode(iclamp,7*i+128,new TGeoCombiTrans(x,y,-38.25,rot12));
+ v9->AddNode(iclamp,7*i+129,new TGeoCombiTrans(x,y,-80.25,rot12));
+ v9->AddNode(iclamp,7*i+130,new TGeoCombiTrans(x,y,-122.25,rot12));
+ v9->AddNode(iclamp,7*i+131,new TGeoCombiTrans(x,y,-164.25,rot12));
+ v9->AddNode(iclamp,7*i+132,new TGeoCombiTrans(x,y,-206.25,rot12));
+ v9->AddNode(iclamp,7*i+133,new TGeoCombiTrans(x,y,-246.75,rot12));
+ //
+ //--------------------------
// outer rods
r=254.25;
x=r * TMath::Cos(angle);
//
// outer rod holder + outer left plug
//
+
TGeoRotation *rot33 = new TGeoRotation();
rot33->RotateZ(-90+i*20.);
//
//
if(i==3){//resistor rod outer
- rotrod.RotateZ(90.+angle);
+ rotrod.RotateZ(90.+i*20.);
rotrod1[1]= new TGeoRotation();
rotpos[1]= new TGeoRotation();
- rotrod1[1]->RotateZ(90.+angle);
+ rotrod1[1]->RotateZ(90.+i*20.);
*rotpos[1] = refl*rotrod;//rotation+reflection
v9->AddNode(tpcrrod,3,new TGeoCombiTrans(x,y, z, rotrod1[1])); //A
v9->AddNode(tpcrrod,4,new TGeoCombiTrans(x,y, -z, rotpos[1])); //C
v9->AddNode(tpcmrod,i+55,new TGeoCombiTrans(x,y,-z,ref));//muon
}
if(i==15){
- v9->AddNode(hvrv,1,new TGeoTranslation(x,y,z+0.7)); //hv->A-side only
-
-
+ v9->AddNode(hvrv,1,new TGeoTranslation(x,y,z+0.7)); //hv->A-side only
}
+ //
+ // outer clamps
+ //
+ r=256.9;
+ x=r * TMath::Cos(angle);
+ y=r * TMath::Sin(angle);
+ rot12= new TGeoRotation();
+ rot12->RotateZ(i*20.);
+ //
+ //A-side
+ v9->AddNode(oclamp,7*i+1,new TGeoCombiTrans(x,y,5.25,rot12));
+ v9->AddNode(oclamp,7*i+2,new TGeoCombiTrans(x,y,38.25,rot12));
+ v9->AddNode(oclamp,7*i+3,new TGeoCombiTrans(x,y,80.25,rot12));
+ v9->AddNode(oclamp,7*i+4,new TGeoCombiTrans(x,y,122.25,rot12));
+ v9->AddNode(oclamp,7*i+5,new TGeoCombiTrans(x,y,164.25,rot12));
+ v9->AddNode(oclamp,7*i+6,new TGeoCombiTrans(x,y,206.25,rot12));
+ v9->AddNode(oclamp,7*i+7,new TGeoCombiTrans(x,y,246.75,rot12));
+ //C-side
+ v9->AddNode(oclamp,7*i+127,new TGeoCombiTrans(x,y,-5.25,rot12));
+ v9->AddNode(oclamp,7*i+128,new TGeoCombiTrans(x,y,-38.25,rot12));
+ v9->AddNode(oclamp,7*i+129,new TGeoCombiTrans(x,y,-80.25,rot12));
+ v9->AddNode(oclamp,7*i+130,new TGeoCombiTrans(x,y,-122.25,rot12));
+ v9->AddNode(oclamp,7*i+131,new TGeoCombiTrans(x,y,-164.25,rot12));
+ v9->AddNode(oclamp,7*i+132,new TGeoCombiTrans(x,y,-206.25,rot12));
+ v9->AddNode(oclamp,7*i+133,new TGeoCombiTrans(x,y,-246.75,rot12));
+
} //end of rods positioning
TGeoVolume *alice = gGeoManager->GetVolume("ALIC");
AliTPC::Init();
- fIdSens=gMC->VolId("TPC_Strip"); // one strip is always selected...
+ fIdSens=TVirtualMC::GetMC()->VolId("TPC_Strip"); // one strip is always selected...
- fIDrift=gMC->VolId("TPC_Drift");
+ fIDrift=TVirtualMC::GetMC()->VolId("TPC_Drift");
fSecOld=-100; // fake number
- gMC->SetMaxNStep(-30000); // max. number of steps increased
+ TVirtualMC::GetMC()->SetMaxNStep(-120000); // max. number of steps increased
if (fPrimaryIonisation) {
// for FLUKA
- gMC->Gstpar(idtmed[2],"PRIMIO_E", 20.77); // 1st ionisation potential
+ TVirtualMC::GetMC()->Gstpar(idtmed[2],"PRIMIO_E", 20.77); // 1st ionisation potential
- gMC->Gstpar(idtmed[2],"PRIMIO_N", 14.35);
- gMC->Gstpar(idtmed[2],"LOSS", 14); // specific energy loss
- gMC->Gstpar(idtmed[2],"STRA",4);
+ TVirtualMC::GetMC()->Gstpar(idtmed[2],"PRIMIO_N", 14.35);
+ TVirtualMC::GetMC()->Gstpar(idtmed[2],"LOSS", 14); // specific energy loss
+ TVirtualMC::GetMC()->Gstpar(idtmed[2],"STRA",4);
}
// specific energy loss for geant3 is now defined in galice.cuts
AliDebug(1,"*** TPC version 2 initialized ***");
- AliDebug(1,Form("Maximum number of steps = %d",gMC->GetMaxNStep()));
+ AliDebug(1,Form("Maximum number of steps = %d",TVirtualMC::GetMC()->GetMaxNStep()));
//
//
// parameters used for the energy loss calculations
//
- const Float_t kprim = 14.35; // number of primary collisions per 1 cm
- const Float_t kpoti = 20.77e-9; // first ionization potential for Ne/CO2
- const Float_t kwIon = 35.97e-9; // energy for the ion-electron pair creation
+ //const Float_t kprim = 14.35; // number of primary collisions per 1 cm
+ //const Float_t kpoti = 20.77e-9; // first ionization potential for Ne/CO2
+ //const Float_t kwIon = 35.97e-9; // energy for the ion-electron pair creation
+ const Float_t kScalewIonG4 = 0.85; // scale factor to tune kwIon for Geant4
+ const Float_t kFanoFactorG4 = 0.7; // parameter for smearing the number of ionizations (nel) using Geant4
const Int_t kMaxDistRef =15; // maximal difference between 2 stored references
+ Float_t prim = fTPCParam->GetNprim();
+ Float_t poti = fTPCParam->GetFpot();
+ Float_t wIon = fTPCParam->GetWmean();
const Float_t kbig = 1.e10;
vol[1]=0; // preset row number to 0
//
- if (!fPrimaryIonisation) gMC->SetMaxStep(kbig);
+ if (!fPrimaryIonisation) TVirtualMC::GetMC()->SetMaxStep(kbig);
- if(!gMC->IsTrackAlive()) return; // particle has disappeared
+ if(!TVirtualMC::GetMC()->IsTrackAlive()) return; // particle has disappeared
- Float_t charge = gMC->TrackCharge();
+ Float_t charge = TVirtualMC::GetMC()->TrackCharge();
if(TMath::Abs(charge)<=0.) return; // take only charged particles
// check the sensitive volume
- id = gMC->CurrentVolID(copy); // vol ID and copy number (starts from 1!)
+ id = TVirtualMC::GetMC()->CurrentVolID(copy); // vol ID and copy number (starts from 1!)
if(id != fIDrift && id != fIdSens) return; // not in the sensitive folume
if ( fPrimaryIonisation && id == fIDrift ) {
- Double_t rnd = gMC->GetRandom()->Rndm();
- gMC->SetMaxStep(0.2+(2.*rnd-1.)*0.05); // 2 mm +- rndm*0.5mm step
+ Double_t rnd = TVirtualMC::GetMC()->GetRandom()->Rndm();
+ TVirtualMC::GetMC()->SetMaxStep(0.2+(2.*rnd-1.)*0.05); // 2 mm +- rndm*0.5mm step
}
- //if ( fPrimaryIonisation && id == fIDrift && gMC->IsTrackEntering()) {
- // gMC->SetMaxStep(0.2); // 2 mm
+ //if ( fPrimaryIonisation && id == fIDrift && TVirtualMC::GetMC()->IsTrackEntering()) {
+ // TVirtualMC::GetMC()->SetMaxStep(0.2); // 2 mm
//}
- gMC->TrackPosition(p);
+ TVirtualMC::GetMC()->TrackPosition(p);
Double_t r = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
//
// track is in the sensitive strip
if(id == fIdSens){
// track is entering the strip
- if (gMC->IsTrackEntering()){
+ if (TVirtualMC::GetMC()->IsTrackEntering()){
Int_t totrows = fTPCParam->GetNRowLow()+fTPCParam->GetNRowUp();
vol[1] = (copy<=totrows) ? copy-1 : copy-1-totrows;
// row numbers are autonomous for lower and upper sectors
// lower sector, row 0, because Jouri wants to have this
- gMC->TrackMomentum(p);
+ TVirtualMC::GetMC()->TrackMomentum(p);
hits[0]=p[0];
hits[1]=p[1];
hits[2]=p[2];
hits[3]=0.; // this hit has no energy loss
// Get also the track time for pileup simulation
- hits[4]=gMC->TrackTime();
+ hits[4]=TVirtualMC::GetMC()->TrackTime();
AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol,hits);
}
//
- gMC->TrackPosition(p);
+ TVirtualMC::GetMC()->TrackPosition(p);
hits[0]=p[0];
hits[1]=p[1];
hits[2]=p[2];
hits[3]=0.; // this hit has no energy loss
// Get also the track time for pileup simulation
- hits[4]=gMC->TrackTime();
+ hits[4]=TVirtualMC::GetMC()->TrackTime();
AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol,hits);
//-----------------------------------------------------------------
// charged particle is in the sensitive drift volume
//-----------------------------------------------------------------
- if(gMC->TrackStep() > 0) {
+ if(TVirtualMC::GetMC()->TrackStep() > 0) {
Int_t nel=0;
if (!fPrimaryIonisation) {
- nel = (Int_t)(((gMC->Edep())-kpoti)/kwIon) + 1;
+ nel = (Int_t)(((TVirtualMC::GetMC()->Edep())-poti)/wIon) + 1;
}
else {
+
+ /*
static Double_t deForNextStep = 0.;
// Geant4 (the meaning of Edep as in Geant3) - wrong
- //nel = (Int_t)(((gMC->Edep())-kpoti)/kwIon) + 1;
+ //nel = (Int_t)(((TVirtualMC::GetMC()->Edep())-poti)/wIon) + 1;
// Geant4 (the meaning of Edep as in Geant3) - NEW
- Double_t eAvailable = gMC->Edep() + deForNextStep;
- nel = (Int_t)(eAvailable/kwIon);
- deForNextStep = eAvailable - nel*kwIon;
+ Double_t eAvailable = TVirtualMC::GetMC()->Edep() + deForNextStep;
+ nel = (Int_t)(eAvailable/wIon);
+ deForNextStep = eAvailable - nel*wIon;
+ */
+
+ //new Geant4-approach
+ Double_t meanIon = TVirtualMC::GetMC()->Edep()/(wIon*kScalewIonG4);
+ nel = (Int_t) ( kFanoFactorG4*AliMathBase::Gamma(meanIon/kFanoFactorG4)); // smear nel using gamma distr w mean = meanIon and variance = meanIon/kFanoFactorG4
}
nel=TMath::Min(nel,300); // 300 electrons corresponds to 10 keV
//
- gMC->TrackPosition(p);
+ TVirtualMC::GetMC()->TrackPosition(p);
hits[0]=p[0];
hits[1]=p[1];
hits[2]=p[2];
// if (fHitType&&2){
if(fHitType){
- gMC->TrackMomentum(p);
+ TVirtualMC::GetMC()->TrackMomentum(p);
Float_t momentum = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
Float_t precision = (momentum>0.1) ? 0.002 :0.01;
fTrackHits->SetHitPrecision(precision);
}
// Get also the track time for pileup simulation
- hits[4]=gMC->TrackTime();
+ hits[4]=TVirtualMC::GetMC()->TrackTime();
AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol,hits);
if (fDebugStreamer){
// function CreateDebugStremer() to be called in the Config.C macro
// if you want to enable it
// By default debug streaemer is OFF
- Float_t edep = gMC->Edep();
- Float_t tstep = gMC->TrackStep();
- Int_t pid=gMC->TrackPid();
+ Float_t edep = TVirtualMC::GetMC()->Edep();
+ Float_t tstep = TVirtualMC::GetMC()->TrackStep();
+ Int_t pid=TVirtualMC::GetMC()->TrackPid();
(*fDebugStreamer)<<"hit"<<
"x="<<hits[0]<< // hit position
"y="<<hits[1]<<
TLorentzVector mom;
// below is valid only for Geant3 (fPromaryIonisation not set)
if(!fPrimaryIonisation){
- gMC->TrackMomentum(mom);
+ TVirtualMC::GetMC()->TrackMomentum(mom);
Float_t ptot=mom.Rho();
- Float_t betaGamma = ptot/gMC->TrackMass();
-
- Int_t pid=gMC->TrackPid();
- if((pid==kElectron || pid==kPositron) && ptot > 0.002)
- {
- pp = kprim*1.58; // electrons above 20 MeV/c are on the plateau!
- }
- else
- {
-
- betaGamma = TMath::Max(betaGamma,(Float_t)7.e-3); // protection against too small bg
- pp=kprim*AliMathBase::BetheBlochAleph(betaGamma);
-
- }
+ Float_t betaGamma = ptot/TVirtualMC::GetMC()->TrackMass();
+
+ //Int_t pid=TVirtualMC::GetMC()->TrackPid();
+ // if((pid==kElectron || pid==kPositron) && ptot > 0.002)
+ // {
+ // pp = prim*1.58; // electrons above 20 MeV/c are on the plateau!
+ // }
+ // else
+ // {
+
+ betaGamma = TMath::Max(betaGamma,(Float_t)7.e-3); // protection against too small bg
+ TVectorD *bbpar = fTPCParam->GetBetheBlochParameters(); //get parametrization from OCDB
+ pp=prim*AliMathBase::BetheBlochAleph(betaGamma,(*bbpar)(0),(*bbpar)(1),(*bbpar)(2),(*bbpar)(3),(*bbpar)(4));
+ // }
- Double_t rnd = gMC->GetRandom()->Rndm();
+ Double_t rnd = TVirtualMC::GetMC()->GetRandom()->Rndm();
- gMC->SetMaxStep(-TMath::Log(rnd)/pp);
+ TVirtualMC::GetMC()->SetMaxStep(-TMath::Log(rnd)/pp);
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff