1 // **************************************************************************
2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 // * Author: The ALICE Off-line Project. *
5 // * Contributors are mentioned in the code where appropriate. *
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 purpose. It is *
13 // * provided "as is" without express or implied warranty. *
14 // **************************************************************************
18 #include <TArrayD.h> //new constructor
20 #include <TGeoManager.h> //CreateGeometry()
21 #include <TGeoVolume.h> //CreateGeometry()
22 #include <TVirtualMC.h> //->gMC in StepManager
23 #include <TPDGCode.h> //StepHistory
24 #include <TClonesArray.h>
25 #include "AliRun.h" //CreateMaterials()
26 #include "AliMC.h" //CreateMaterials()
27 #include "AliMagF.h" //CreateMaterials()
28 #include "AliCDBManager.h" //CreateMaterials()
29 #include "AliCDBEntry.h" //CreateMaterials()
30 #include "AliITSupgrade.h" //class header
31 #include "AliITShit.h" // StepManager()
32 #include "AliITSDigitUpgrade.h"
33 #include "AliTrackReference.h" // StepManager()
34 #include "AliITSDetTypeSim.h"
35 ClassImp(AliITSupgrade)
36 //__________________________________________________________________________________________________
37 AliITSupgrade::AliITSupgrade():
58 //__________________________________________________________________________________________________
59 AliITSupgrade::AliITSupgrade(const char *name,const char *title, TArrayD widths, TArrayD radii,TArrayD halfLengths, TArrayD radiiCu, TArrayD widthsCu, TArrayS copper,Bool_t bp,Double_t radiusBP, Double_t widthBP, Double_t halfLengthBP):
69 fHalfLengthBP(halfLengthBP),
70 fNlayers(widths.GetSize()),
71 fHalfLength(halfLengths),
80 for(Int_t i=0;i<fNlayers;i++){
82 fWidths.AddAt(widths.At(i),i);
84 fRadii.AddAt(radii.At(i),i);
85 AliDebug(1,"Creating Volume");
91 AliITSupgrade::~AliITSupgrade(){
93 if(fSdigits) {fSdigits->Delete(); delete fSdigits;}
94 if(fDigits) {fDigits->Delete(); delete fDigits;}
99 //_________________________________________________________________________________________________
100 void AliITSupgrade::AddAlignableVolumes()const
106 //__________________________________________________________________________________________________
108 void AliITSupgrade::CreateMaterials()
111 // Definition of ITS materials
114 AliInfo("Start ITS materials");
115 //data from PDG booklet 2002 density [gr/cm^3] rad len [cm] abs len [cm]
116 Float_t aAir[4]={12,14,16,36} , zAir[4]={6,7,8,18} , wAir[4]={0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; Int_t nAir=4;//mixture
117 Float_t aBe = 9.012 , zBe = 4 , dBe = 1.848 , radBe = 65.19/dBe , absBe = 77.8/dBe ; // UPGRADE -> beryllium beampipe
118 Float_t aSi = 28.085 , zSi = 14 , dSi = 2.329 , radSi = 21.82/dSi , absSi = 108.4/dSi ; // UPGRADE -> Si tube
119 Float_t aCu = 63.54, zCu = 29 , dCu= 8.96 , radCu = 12.86/dCu, absCu= 137.3/dCu;//Upgrade -> Copper Tube
121 Int_t matId=0; //tmp material id number
122 Int_t unsens = 0, sens=1; //sensitive or unsensitive medium
123 Int_t itgfld = 3; //type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
124 Float_t maxfld = 5.; //max field value
125 Float_t tmaxfd = -10.0; //max deflection angle due to magnetic field in one step
126 Float_t deemax = - 0.2; //max fractional energy loss in one step
127 Float_t stemax = - 0.1; //max step allowed [cm]
128 Float_t epsil = 0.001; //abs tracking precision [cm]
129 Float_t stmin = - 0.001; //min step size [cm] in continius process transport, negative value: choose it automatically
130 Float_t tmaxfdSi = 0.1; // .10000E+01; // Degree
131 Float_t stemaxSi = 0.0075; // .10000E+01; // cm
132 Float_t deemaxSi = 0.1; // 0.30000E-02; // Fraction of particle's energy 0<deemax<=1
133 Float_t epsilSi = 1.0E-4;// .10000E+01;
134 Float_t stminSi = 0.0; // cm "Default value used"
136 Float_t epsilBe = .001; // Tracking precision,
137 Float_t stemaxBe = -0.01; // Maximum displacement for multiple scat
138 Float_t tmaxfdBe = -20.; // Maximum angle due to field deflection
139 Float_t deemaxBe = -.3; // Maximum fractional energy loss, DLS
140 Float_t stminBe = -.8;
143 AliMixture(++matId,"UpgradeAir" ,aAir ,zAir ,dAir ,nAir ,wAir );
144 AliMedium(kAir ,"UpgradeAir" ,matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
146 AliMaterial(++matId,"UpgradeSi" ,aSi ,zSi ,dSi ,radSi ,absSi );
147 AliMedium(kSi ,"UpgradeSi" , matId, sens, itgfld, maxfld, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
149 AliMaterial(++matId,"UpgradeBe" ,aBe ,zBe ,dBe ,radBe ,absBe );
150 AliMedium(kBe ,"UpgradeBe" , matId, unsens, itgfld, maxfld, tmaxfdBe, stemaxBe, deemaxBe, epsilBe, stminBe);
152 AliMaterial(++matId, "UpgradeCu", aCu, zCu, dCu, radCu, absCu);
153 AliMedium(kCu, "UpgradeCu", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
155 AliInfo("End ITS materials");
157 }//void AliITS::CreateMaterials()
158 //__________________________________________________________________________________________________
159 void AliITSupgrade::CreateGeometry()
162 //Creates detailed geometry simulation (currently GEANT volumes tree)
164 AliInfo("Start ITS upgrade preliminary version building");
165 if(!gMC->IsRootGeometrySupported()) return;
166 TGeoVolumeAssembly *vol= CreateVol();
167 gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
168 AliInfo("Stop ITS upgrade preliminary version building");
170 //__________________________________________________________________________________________________
171 void AliITSupgrade::Init()
173 // This method defines ID for sensitive volumes, i.e. such geometry volumes for which there are if(gMC->CurrentVolID()==XXX)
174 // statements in StepManager()
177 AliDebug(1,"Init ITS upgrade preliminary version.");
179 //__________________________________________________________________________________________________
180 void AliITSupgrade::StepManager()
182 // Full Step Manager.
185 // StepHistory(); return; //uncomment to print tracks history
186 // StepCount(); return; //uncomment to count photons
187 if(!(this->IsActive())) return;
188 if(!(gMC->TrackCharge())) return;
189 TString volumeName=gMC->CurrentVolName();
190 if(gMC->IsTrackExiting() && !volumeName.Contains("Cu") && !volumeName.Contains("Be")) {
191 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kITS);
192 } // if Outer ITS mother Volume
194 static TLorentzVector position, momentum; // Saves on calls to construtors
195 static AliITShit hit;// Saves on calls to constructors
202 if(gMC->IsTrackInside()) status += 1;
203 if(gMC->IsTrackEntering()) status += 2;
204 if(gMC->IsTrackExiting()) status += 4;
205 if(gMC->IsTrackOut()) status += 8;
206 if(gMC->IsTrackDisappeared()) status += 16;
207 if(gMC->IsTrackStop()) status += 32;
208 if(gMC->IsTrackAlive()) status += 64;
211 // Fill hit structure.
213 TString volname = gMC->CurrentVolName();
214 if(volname.Contains("Cu"))return;
215 if(volname.Contains("Be"))return;
216 //AliInfo(Form("volume name %s ",volname.Data()));
217 volname.Remove(0,12); // remove letters to get the layer number
218 //AliInfo(Form("volume name %s ",volname.Data()));
219 hit.SetModule(volname.Atoi()); // this will be the layer, not the module
220 hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
222 gMC->TrackPosition(position);
223 gMC->TrackMomentum(momentum);
224 hit.SetPosition(position);
226 hit.SetTime(gMC->TrackTime());
227 hit.SetMomentum(momentum);
228 hit.SetStatus(status);
229 hit.SetEdep(gMC->Edep());
230 hit.SetShunt(GetIshunt());
231 if(gMC->IsTrackEntering()){
232 hit.SetStartPosition(position);
233 hit.SetStartTime(gMC->TrackTime());
234 hit.SetStartStatus(status);
235 return; // don't save entering hit.
237 // Fill hit structure with this new hit.
239 new((*fHits)[fNhits++]) AliITShit(hit); // Use Copy Construtor.
240 // Save old position... for next hit.
241 hit.SetStartPosition(position);
242 hit.SetStartTime(gMC->TrackTime());
243 hit.SetStartStatus(status);
247 //__________________________________________________________________________________________________
248 TGeoVolumeAssembly * AliITSupgrade::CreateVol()
250 TGeoVolumeAssembly *vol = new TGeoVolumeAssembly("ITSupgrade");
251 TGeoMedium *si =gGeoManager->GetMedium("ITS_UpgradeSi");
252 TGeoMedium *cu =gGeoManager->GetMedium("ITS_UpgradeCu");
253 TGeoMedium *be =gGeoManager->GetMedium("ITS_UpgradeBe");
254 for(Int_t ivol=0;ivol<fNlayers;ivol++){
255 TGeoVolume *layer=gGeoManager->MakeTube(Form("LayerSilicon%i",ivol),si , fRadii.At(ivol) , fRadii.At(ivol)+fWidths.At(ivol) , fHalfLength.At(ivol)); //upgraded situation
257 TGeoVolume *layerCu=gGeoManager->MakeTube(Form("LayerCu%i",ivol),cu , fRadiiCu.At(ivol) , fRadiiCu.At(ivol)+fWidthsCu.At(ivol) , fHalfLength.At(ivol) ); //upgraded situation
260 vol ->AddNode(layer,ivol);
261 if(fCopper.At(ivol)){
262 vol->AddNode(layerCu,ivol);
265 TGeoVolume *beampipe=gGeoManager->MakeTube("BeamPipe", be , fRadiusBP , fRadiusBP+ fWidthBP , fHalfLengthBP ); //upgraded situation
267 if(fBeampipe) vol->AddNode(beampipe,0);
271 //_________________________________________________________________________________________________
272 void AliITSupgrade::SetFullSegmentation(TArrayD xsize,TArrayD zsize){
273 TFile *file= TFile::Open("Segmentation.root","recreate");
274 file->WriteObjectAny(&xsize,"TArrayD","CellSizeX");
275 file->WriteObjectAny(&zsize,"TArrayD","CellSizeZ");
278 //_________________________________________________________________________________________________
279 void AliITSupgrade::StepHistory()
281 // This methode is invoked from StepManager() in order to print out
283 const char *sParticle;
284 switch(gMC->TrackPid()){
285 case kProton: sParticle="PROTON" ;break;
286 case kNeutron: sParticle="neutron" ;break;
287 case kGamma: sParticle="gamma" ;break;
288 case kPi0: sParticle="Pi0" ;break;
289 case kPiPlus: sParticle="Pi+" ;break;
290 case kPiMinus: sParticle="Pi-" ;break;
291 case kElectron: sParticle="electron" ;break;
292 default: sParticle="not known" ;break;
295 TString flag="fanny combination";
296 if(gMC->IsTrackAlive())
297 if(gMC->IsTrackEntering()) flag="enters to";
298 else if(gMC->IsTrackExiting()) flag="exits from";
299 else if(gMC->IsTrackInside()) flag="inside";
300 else if(gMC->IsTrackStop()) flag="stopped in";
303 TString path=gMC->CurrentVolName(); path.Prepend("-");path.Prepend(gMC->CurrentVolOffName(1));//current volume and his mother are always there
304 vid=gMC->CurrentVolOffID(2,copy); if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
305 vid=gMC->CurrentVolOffID(3,copy); if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
307 AliDebug(10, Form("Step %i: %s (%i) %s %s m=%.6f GeV q=%.1f dEdX=%.4f Etot=%.4f",iStepN,sParticle,gMC->TrackPid(),flag.Data(),path.Data(),gMC->TrackMass(),gMC->TrackCharge(),gMC->Edep()*1e9,gMC->Etot()));
309 Double_t gMcTrackPos[3]; gMC->TrackPosition(gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2]);
310 Double_t gMcTrackPosLoc[3]; gMC->Gmtod(gMcTrackPos,gMcTrackPosLoc,1);
311 AliDebug(10,Form("gMC Track Position (MARS) x: %5.3lf, y: %5.3lf, z: %5.3lf (r: %5.3lf) ---> (LOC) x: %5.3f, y: %5.3f, z: %5.3f",gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2],TMath::Sqrt(gMcTrackPos[0]*gMcTrackPos[0]+gMcTrackPos[1]*gMcTrackPos[1]+gMcTrackPos[2]*gMcTrackPos[2]),gMcTrackPosLoc[0],gMcTrackPosLoc[1],gMcTrackPosLoc[2]));
314 AliDebug(10,Form("Step %i: tid=%i flags alive=%i disap=%i enter=%i exit=%i inside=%i out=%i stop=%i new=%i",
315 iStepN, gAlice->GetMCApp()->GetCurrentTrackNumber(),
316 gMC->IsTrackAlive(), gMC->IsTrackDisappeared(),gMC->IsTrackEntering(), gMC->IsTrackExiting(),
317 gMC->IsTrackInside(),gMC->IsTrackOut(), gMC->IsTrackStop(), gMC->IsNewTrack()));
319 Float_t a,z,den,rad,abs; a=z=den=rad=abs=-1;
320 Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs);
321 AliDebug(10, Form("Step %i: mid=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f\n\n",iStepN,mid,a,z,den,rad,abs));
323 TArrayI proc; gMC->StepProcesses(proc);
324 AliDebug(1,"Processes in this step:");
325 for ( int i = 0 ; i < proc.GetSize(); i++)
327 AliDebug(10, Form("%s",TMCProcessName[proc.At(i)]));
329 AliDebug(1,"End process list");
332 //______________________________________________________
333 void AliITSupgrade::Hits2SDigits(){
335 // Interface method ivoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
338 AliDebug(1,"Start Hits2SDigits.");
340 if(!fSegmentation)fSegmentation=new AliITSsegmentationUpgrade();
342 if(!GetLoader()->TreeH()) {GetLoader()->LoadHits();}
343 if(!GetLoader()->TreeS())
345 for(Int_t iEvt=0;iEvt < GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvt++){//events loop
346 GetLoader()->GetRunLoader()->GetEvent(iEvt); //get next event
348 GetLoader()->MakeTree("S");
352 Int_t nSdigit[10]; for(Int_t i=0;i<10;i++) nSdigit[i] =0;
353 for(Int_t iEnt=0;iEnt<GetLoader()->TreeH()->GetEntries();iEnt++){//prims loop
354 GetLoader()->TreeH()->GetEntry(iEnt);
355 Hit2SumDig(Hits(),SDigitsList(),nSdigit);//convert this hit to list of sdigits
357 GetLoader()->TreeS()->Fill();
358 GetLoader()->WriteSDigits("OVERWRITE");
361 GetLoader()->UnloadHits();
362 GetLoader()->UnloadSDigits();
363 AliDebug(1,"Stop Hits2SDigits.");
366 //____________________________
367 void AliITSupgrade::Hit2SumDig(TClonesArray *hits,TObjArray *pSDig, Int_t *nSdigit)
369 // Adds sdigits of this hit to the list
372 AliDebug(1,"Start Hit2SumDig");
376 AliDebug(1,"Segmentation Not inizialized!!");
380 TClonesArray *pSdigList[10];
382 for(Int_t i=0;i<fNlayers;i++){
383 pSdigList[i]=(TClonesArray*)(*pSDig)[i];
384 if(pSdigList[i]->GetEntries()!=0) AliErrorClass("Some of sdigits lists is not empty"); //in principle those lists should be empty
387 AliInfo(Form("Number of hits %i ",hits->GetEntries()));
388 for(Int_t iHit=0;iHit<hits->GetEntries();iHit++){ //hits loop
389 AliITShit *hit = (AliITShit*)hits->At(iHit);
391 if(!fSegmentation->GlobalToDet(hit->GetModule(),hit->GetXG(),hit->GetYG(),hit->GetZG(),xz[0],xz[1])) continue;
392 AliITSDigitUpgrade digit;
393 digit.SetSignal(hit->GetIonization());
394 digit.SetNelectrons(hit->GetIonization()/(3.62*1e-09));
395 digit.SetLayer(hit->GetModule());
396 digit.SetTrackID(hit->GetTrack());
399 if(fSegmentation->GetCellSizeX(hit->GetModule())!=0) xpix =(Int_t) (xz[0]/ fSegmentation->GetCellSizeX(hit->GetModule()));
401 Int_t zpix = 999; // shift at the next line to have zpix always positive. Such numbers are used to build the Pixel Id in the layer (> 0!)
402 if(fSegmentation->GetCellSizeZ(hit->GetModule())!=0){
403 zpix =(Int_t)((xz[1]+fHalfLength.At(hit->GetModule())) / fSegmentation->GetCellSizeZ(hit->GetModule()));
405 digit.SetPixId(xpix,zpix);
408 new((*pSdigList[hit->GetModule()])[nSdigit[hit->GetModule()]++]) AliITSDigitUpgrade(digit);
411 AliDebug(1,"Stop Hit2SumDig.");
413 //_______________________________________________________________________________________________
414 void AliITSupgrade::MakeBranch(Option_t *option){
415 //Create Tree branches
416 AliDebug(1,Form("Start with option= %s.",option));
418 const Int_t kBufSize = 4000;
420 const char *cH = strstr(option,"H");
421 const char *cD = strstr(option,"D");
422 const char *cR = strstr(option,"R");
423 const char *cS = strstr(option,"S");
425 if(cH&&fLoader->TreeH()){
428 MakeBranchInTree(fLoader->TreeH(),"ITSupgrade",&fHits,kBufSize,0);
432 if(cS&&fLoader->TreeS()){
434 for(Int_t i=0;i<fNlayers;i++) MakeBranchInTree(fLoader->TreeS(),Form("Layer%d",i),&((*fSdigits)[i]),kBufSize,0);
438 if(cD&&fLoader->TreeD()){
440 for(Int_t i=0;i<fNlayers;i++) MakeBranchInTree(fLoader->TreeD(),Form("Layer%d",i),&((*fDigits)[i]),kBufSize,0);
443 if(cR&&fLoader->TreeR()){
445 for(Int_t i=0;i<fNlayers;i++) MakeBranchInTree(fLoader->TreeR(),Form("Layer%d",i),&((*fClusters)[i]),kBufSize,0);
450 //____________________________________________________________________________________________________
451 void AliITSupgrade::SetTreeAddress()
453 //Set branch address for the Hits and Digits Tree.
454 AliDebug(1,"Start.");
455 if(fLoader->TreeH() && fLoader->TreeH()->GetBranch("ITSupgrade" )){
457 fLoader->TreeH()->SetBranchAddress("ITSupgrade",&fHits);
463 if(fLoader->TreeS() && fLoader->TreeS()->GetBranch("Layer0" )){
465 for(int i=0;i<fNlayers;i++){
466 fLoader->TreeS()->SetBranchAddress(Form("Layer%d",i),&((*fSdigits)[i]));
471 if(fLoader->TreeD() && fLoader->TreeD()->GetBranch("Layer0")){
473 for(int i=0;i<fNlayers;i++){
474 fLoader->TreeD()->SetBranchAddress(Form("Layer%d",i),&((*fDigits)[i]));
477 if(fLoader->TreeR() && fLoader->TreeR()->GetBranch("Layer0")){
479 for(int i=0;i<fNlayers;i++){
480 fLoader->TreeR()->SetBranchAddress(Form("Layer%d",i),&((*fClusters)[i]));