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 Revision 1.19 2002/09/09 17:23:28 nilsen
19 Minor changes in support of changes to AliITSdigitS?D class'.
21 Revision 1.18 2002/08/21 22:11:13 nilsen
22 Debug output now settable via a DEBUG flag.
24 Revision 1.17 2002/07/16 17:00:17 barbera
25 Fixes added to make the slow simulation running with the current HEAD (from M. Masera)
27 Revision 1.16 2002/06/19 16:02:22 hristov
28 Division by zero corrected
30 Revision 1.15 2002/03/15 17:32:14 nilsen
31 Reintroduced SDigitization, and Digitization from SDigits, along with
32 functions InitSimulationModule, and FinishSDigitizModule.
34 Revision 1.14 2001/11/23 13:04:07 barbera
35 Some protection added in case of high multiplicity
37 Revision 1.13 2001/11/13 11:13:24 barbera
38 A protection against tracks with the same entrance and exit has been made more strict
40 Revision 1.12 2001/10/04 22:44:31 nilsen
41 Major changes in supppor of PreDigits (SDigits). Changes made with will make
42 it easier to suppor expected changes in AliITSHit class. Added use of new
43 class AliITSpList. Both SPD and SDD have added effects of Dead Channels. Both
44 of these will require addtional work as data bases of detectors and the like
53 #include <TParticle.h>
57 #include "AliITShit.h"
58 #include "AliITSdigit.h"
59 #include "AliITSmodule.h"
60 #include "AliITSMapA2.h"
61 #include "AliITSpList.h"
62 #include "AliITSsimulationSPD.h"
63 #include "AliITSsegmentation.h"
64 #include "AliITSresponse.h"
65 #include "AliITSsegmentationSPD.h"
66 #include "AliITSresponseSPD.h"
70 ClassImp(AliITSsimulationSPD)
71 ////////////////////////////////////////////////////////////////////////
73 // Written by Rocco Caliandro
74 // from a model developed with T. Virgili and R.A. Fini
77 // AliITSsimulationSPD is the simulation of SPDs
79 //______________________________________________________________________
80 AliITSsimulationSPD::AliITSsimulationSPD(){
81 // Default constructor
94 //______________________________________________________________________
95 AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg,
96 AliITSresponse *resp) {
97 // Standard constructor
109 Init((AliITSsegmentationSPD*)seg,(AliITSresponseSPD*)resp);
111 //______________________________________________________________________
112 void AliITSsimulationSPD::Init(AliITSsegmentationSPD *seg,
113 AliITSresponseSPD *resp) {
114 // Initilizes the variables of AliITSsimulation SPD.
119 fMapA2 = new AliITSMapA2(fSegmentation);
120 fpList = new AliITSpList(GetNPixelsZ()+1,GetNPixelsX()+1);
122 fResponse->Thresholds(fThresh,fSigma);
123 fResponse->GetNoiseParam(fCouplCol,fCouplRow);
124 fNPixelsZ = fSegmentation->Npz();
125 fNPixelsX = fSegmentation->Npx();
128 //______________________________________________________________________
129 AliITSsimulationSPD::~AliITSsimulationSPD() {
140 //______________________________________________________________________
141 AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
144 if(&source == this) return;
146 this->fMapA2 = source.fMapA2;
147 this->fHis = source.fHis;
149 this->fThresh = source.fThresh;
150 this->fSigma = source.fSigma;
151 this->fCouplCol = source.fCouplCol;
152 this->fCouplRow = source.fCouplRow;
153 this->fNPixelsX = source.fNPixelsX;
154 this->fNPixelsZ = source.fNPixelsZ;
158 //______________________________________________________________________
159 AliITSsimulationSPD& AliITSsimulationSPD::operator=(const AliITSsimulationSPD
161 // Assignment operator
163 if(&source == this) return *this;
165 this->fMapA2 = source.fMapA2;
166 this->fHis = source.fHis;
168 this->fThresh = source.fThresh;
169 this->fSigma = source.fSigma;
170 this->fCouplCol = source.fCouplCol;
171 this->fCouplRow = source.fCouplRow;
172 this->fNPixelsX = source.fNPixelsX;
173 this->fNPixelsZ = source.fNPixelsZ;
177 //______________________________________________________________________
178 void AliITSsimulationSPD::InitSimulationModule(Int_t module,Int_t event){
179 // Creates maps to build the list of tracks for each sumable digit
181 // Int_t module // Module number to be simulated
182 // Int_t event // Event number to be simulated
193 //______________________________________________________________________
194 void AliITSsimulationSPD::FinishSDigitiseModule(){
195 // Does the Sdigits to Digits work
203 SDigitsToDigits(fModule,fpList);
205 //______________________________________________________________________
206 void AliITSsimulationSPD::SDigitiseModule(AliITSmodule *mod, Int_t dummy0,
208 // Sum digitize module
209 if (!(mod->GetNhits())) return; // if module has no hits then no Sdigits.
210 Int_t number = 10000;
211 Int_t *frowpixel = new Int_t[number];
212 Int_t *fcolpixel = new Int_t[number];
213 Double_t *fenepixel = new Double_t[number];
215 fModule = mod->GetIndex();
217 // Array of pointers to store the track index of the digits
218 // leave +1, otherwise pList crashes when col=256, row=192
220 HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,fpList);
222 WriteSDigits(fpList);
231 //______________________________________________________________________
232 void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t dummy0,
234 // digitize module. Also need to digitize modules with only noise.
236 Int_t number = 10000;
237 Int_t *frowpixel = new Int_t[number];
238 Int_t *fcolpixel = new Int_t[number];
239 Double_t *fenepixel = new Double_t[number];
241 // Array of pointers to store the track index of the digits
242 // leave +1, otherwise pList crashes when col=256, row=192
243 fModule = mod->GetIndex();
245 SetFluctuations(fpList,fModule);
247 HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,fpList);
249 // apply mask to SPD module
252 CreateDigit(fModule,fpList);
261 //______________________________________________________________________
262 void AliITSsimulationSPD::SDigitsToDigits(Int_t module,AliITSpList *pList) {
263 // sum digits to Digits.
266 cout << "Entering AliITSsimulatinSPD::SDigitsToDigits for module=";
267 cout << module << endl;
272 SetFluctuations(pList,module);
274 fMapA2->ClearMap(); // since noise is in pList aready. Zero Map so that
275 // noise is not doubled when calling FillMapFrompList.
277 FillMapFrompList(pList);
279 // apply mask to SPD module
282 CreateDigit(module,pList);
287 //______________________________________________________________________
288 void AliITSsimulationSPD::UpdateMapSignal(Int_t row,Int_t col,Int_t trk,
289 Int_t hit,Int_t mod,Double_t ene,
290 AliITSpList *pList) {
291 // updates the Map of signal, adding the energy (ene) released by
294 fMapA2->AddSignal(row,col,ene);
295 pList->AddSignal(row,col,trk,hit,mod,ene);
297 //______________________________________________________________________
298 void AliITSsimulationSPD::UpdateMapNoise(Int_t row,Int_t col,Int_t mod,
299 Double_t ene,AliITSpList *pList) {
300 // updates the Map of noise, adding the energy (ene) give my noise
302 fMapA2->AddSignal(row,col,ene);
303 pList->AddNoise(row,col,mod,ene);
305 //______________________________________________________________________
306 void AliITSsimulationSPD::HitsToAnalogDigits(AliITSmodule *mod,
307 Int_t *frowpixel,Int_t *fcolpixel,
309 AliITSpList *pList) {
310 // Loops over all hits to produce Analog/floting point digits. This
311 // is also the first task in producing standard digits.
313 // loop over hits in the module
314 Int_t hitpos,nhits = mod->GetNhits();
315 for (hitpos=0;hitpos<nhits;hitpos++) {
316 HitToDigit(mod,hitpos,frowpixel,fcolpixel,fenepixel,pList);
317 }// end loop over digits
319 //______________________________________________________________________
320 void AliITSsimulationSPD::HitToDigit(AliITSmodule *mod,Int_t hitpos,
321 Int_t *frowpixel,Int_t *fcolpixel,
322 Double_t *fenepixel,AliITSpList *pList) {
323 // Steering function to determine the digits associated to a given
325 // The digits are created by charge sharing (ChargeSharing) and by
326 // capacitive coupling (SetCoupling). At all the created digits is
327 // associated the track number of the hit (ntrack)
328 Double_t x1l=0.0,y1l=0.0,z1l=0.0,x2l=0.0,y2l=0.0,z2l=0.0;
329 Int_t r1,r2,c1,c2,row,col,npixel = 0;
331 Double_t ene=0.0,etot=0.0;
332 const Float_t kconv = 10000.; // cm -> microns
333 const Float_t kconv1= 0.277e9; // GeV -> electrons equivalent
335 if(!(mod->LineSegmentL(hitpos,x1l,x2l,y1l,y2l,z1l,z2l,etot,ntrack)))return;
337 x2l += x1l; y2l += y1l; z2l += z1l; // Convert to ending coordinate.
338 // positions shifted and converted in microns
339 x1l = x1l*kconv + fSegmentation->Dx()/2.;
340 z1l = z1l*kconv + fSegmentation->Dz()/2.;
341 // positions shifted and converted in microns
342 x2l = x2l*kconv + fSegmentation->Dx()/2.;
343 z2l = z2l*kconv + fSegmentation->Dz()/2.;
344 etot *= kconv1; // convert from GeV to electrons equivalent.
345 Int_t module = mod->GetIndex();
347 // to account for the effective sensitive area
348 // introduced in geometry
349 if (z1l<0 || z1l>fSegmentation->Dz()) return;
350 if (z2l<0 || z2l>fSegmentation->Dz()) return;
351 if (x1l<0 || x1l>fSegmentation->Dx()) return;
352 if (x2l<0 || x2l>fSegmentation->Dx()) return;
354 //Get the col and row number starting from 1
355 // the x direction is not inverted for the second layer!!!
356 fSegmentation->GetPadIxz(x1l, z1l, c1, r1);
357 fSegmentation->GetPadIxz(x2l, z2l, c2, r2);
359 // to account for unexpected equal entrance and
361 if (x1l==x2l) x2l=x2l+x2l*0.1;
362 if (z1l==z2l) z2l=z2l+z2l*0.1;
364 if ((r1==r2) && (c1==c2)){
367 frowpixel[npixel-1] = r1;
368 fcolpixel[npixel-1] = c1;
369 fenepixel[npixel-1] = etot;
372 ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
373 npixel,frowpixel,fcolpixel,fenepixel);
374 } // end if r1==r2 && c1==c2.
376 for (Int_t npix=0;npix<npixel;npix++){
377 row = frowpixel[npix];
378 col = fcolpixel[npix];
379 ene = fenepixel[npix];
380 UpdateMapSignal(row,col,ntrack,hitpos,module,ene,pList);
381 // Starting capacitive coupling effect
382 SetCoupling(row,col,ntrack,hitpos,module,pList);
385 //______________________________________________________________________
386 void AliITSsimulationSPD::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
387 Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
388 Int_t r2,Float_t etot,
389 Int_t &npixel,Int_t *frowpixel,
390 Int_t *fcolpixel,Double_t *fenepixel){
391 // Take into account the geometrical charge sharing when the track
392 // crosses more than one pixel.
396 <img src="picts/ITS/barimodel_2.gif">
399 <font size=+2 color=red>
400 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
406 Float_t xa,za,xb,zb,dx,dz,dtot,refr,refm,refc;
408 Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
409 Int_t dirx,dirz,rb,cb;
410 Int_t flag,flagrow,flagcol;
420 dtot = TMath::Sqrt((dx*dx)+(dz*dz));
421 if (dtot==0.0) dtot = 0.01;
422 dirx = (Int_t) TMath::Sign((Float_t)1,dx);
423 dirz = (Int_t) TMath::Sign((Float_t)1,dz);
425 // calculate the x coordinate of the pixel in the next column
426 // and the z coordinate of the pixel in the next row
429 fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos);
431 Float_t xsize = fSegmentation->Dpx(0);
432 Float_t zsize = fSegmentation->Dpz(r1-1);
434 if (dirx == 1) refr = xpos+xsize/2.;
435 else refr = xpos-xsize/2.;
437 if (dirz == 1) refn = zpos+zsize/2.;
438 else refn = zpos-zsize/2.;
444 // calculate the x coordinate of the intersection with the pixel
445 // in the next cell in row direction
447 refm = dx*((refn - z1l)/dz) + x1l;
449 refm = refr+dirx*xsize;
451 // calculate the z coordinate of the intersection with the pixel
452 // in the next cell in column direction
454 refc = dz*((refr - x1l)/dx) + z1l;
456 refc = refn+dirz*zsize;
463 if ((arefm < arefr) && (arefn < arefc)){
464 // the track goes in the pixel in the next cell in row direction
471 if (rb == r2) flagrow=1;
476 // shift to the pixel in the next cell in row direction
477 Float_t zsizeNext = fSegmentation->Dpz(rb-1);
478 //to account for cell at the borders of the detector
479 if(zsizeNext==0) zsizeNext = zsize;
480 refn += zsizeNext*dirz;
482 // the track goes in the pixel in the next cell in column direction
489 if (cb == c2) flagcol=1;
493 } // end ifaxb > ax2l
495 // shift to the pixel in the next cell in column direction
496 Float_t xsizeNext = fSegmentation->Dpx(cb-1);
497 //to account for cell at the borders of the detector
498 if(xsizeNext==0) xsizeNext = xsize;
499 refr += xsizeNext*dirx;
500 } // end if (arefm < arefr) && (arefn < arefc)
502 //calculate the energy lost in the crossed pixel
503 epar = TMath::Sqrt((xb-xa)*(xb-xa)+(zb-za)*(zb-za));
504 epar = etot*(epar/dtot);
506 //store row, column and energy lost in the crossed pixel
507 frowpixel[npixel] = r1;
508 fcolpixel[npixel] = c1;
509 fenepixel[npixel] = epar;
512 // the exit point of the track is reached
513 if (epar == 0) flag = 1;
514 if ((r1 == r2) && (c1 == c2)) flag = 1;
523 //______________________________________________________________________
524 void AliITSsimulationSPD::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
525 Int_t idhit,Int_t module,
526 AliITSpList *pList) {
527 // Take into account the coupling between adiacent pixels.
528 // The parameters probcol and probrow are the fractions of the
529 // signal in one pixel shared in the two adjacent pixels along
530 // the column and row direction, respectively.
534 <img src="picts/ITS/barimodel_3.gif">
537 <font size=+2 color=red>
538 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
544 Double_t pulse1,pulse2;
545 Float_t couplR=0.0,couplC=0.0;
547 GetCouplings(couplR,couplC);
550 pulse1 = fMapA2->GetSignal(row,col);
552 for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
556 if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
557 pulse1 = fMapA2->GetSignal(row,col);
561 UpdateMapSignal(j1,col,ntrack,idhit,module,pulse1,pList);
565 // loop in column direction
569 if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
570 pulse2 = fMapA2->GetSignal(row,col);
574 UpdateMapSignal(row,j2,ntrack,idhit,module,pulse2,pList);
580 //______________________________________________________________________
581 void AliITSsimulationSPD::CreateDigit(Int_t module,AliITSpList *pList) {
582 // The pixels are fired if the energy deposited inside them is above
583 // the threshold parameter ethr. Fired pixed are interpreted as digits
584 // and stored in the file digitfilename. One also needs to write out
585 // cases when there is only noise (nhits==0).
587 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
589 Int_t size = AliITSdigitSPD::GetNTracks();
590 Int_t * digits = new Int_t[size];
591 Int_t * tracks = new Int_t[size];
592 Int_t * hits = new Int_t[size];
593 Float_t * charges = new Float_t[size];
596 for (Int_t r=1;r<=GetNPixelsZ();r++) {
597 for (Int_t c=1;c<=GetNPixelsX();c++) {
598 // check if the deposited energy in a pixel is above the
600 Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
601 if ( signal > GetThreshold()) {
602 digits[0] = r-1; // digits starts from 0
603 digits[1] = c-1; // digits starts from 0
605 digits[2] = (Int_t) signal; // the signal is stored in
607 for(j1=0;j1<size;j1++){
608 if(j1<pList->GetNEnteries()){
609 tracks[j1] = pList->GetTrack(r,c,j1);
610 hits[j1] = pList->GetHit(r,c,j1);
618 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
620 cout << " CreateSPDDigit mod=" << fModule << " r,c=" << r;
621 cout <<","<<c<< " sig=" << fpList->GetSignalOnly(r,c);
622 cout << " noise=" << fpList->GetNoise(r,c);
623 cout << " Msig="<< signal << " Thres=" << GetThreshold()<<endl;
625 } // end if of threshold condition
634 //______________________________________________________________________
635 void AliITSsimulationSPD::SetFluctuations(AliITSpList *pList,Int_t module) {
636 // Set the electronic noise and threshold non-uniformities to all the
637 // pixels in a detector.
638 // The parameter fSigma is the squared sum of the sigma due to noise
639 // and the sigma of the threshold distribution among pixels.
643 <img src="picts/ITS/barimodel_1.gif">
646 <font size=+2 color=red>
647 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
652 Float_t thr=0.0,sigm=0.0;
653 Double_t signal,sigma;
656 GetThresholds(thr,sigm);
657 sigma = (Double_t) sigm;
658 for(iz=1;iz<=GetNPixelsZ();iz++){
659 for(ix=1;ix<=GetNPixelsX();ix++){
660 signal = sigma*gRandom->Gaus();
661 fMapA2->SetHit(iz,ix,signal);
662 // insert in the label-signal-hit list the pixels fired
664 pList->AddNoise(iz,ix,module,signal);
665 } // end of loop on pixels
666 } // end of loop on pixels
668 //______________________________________________________________________
669 void AliITSsimulationSPD::SetMask() {
670 // Apply a mask to the SPD module. 1% of the pixel channels are
671 // masked. When the database will be ready, the masked pixels
672 // should be read from it.
676 Float_t perc = ((AliITSresponseSPD*)fResponse)->GetFractionDead();
677 // in this way we get the same set of random numbers for all runs.
678 // This is a cluge for now.
679 static TRandom *rnd = new TRandom();
681 totMask= perc*GetNPixelsZ()*GetNPixelsX();
682 for(im=1;im<totMask;im++){
684 ix=(Int_t)(rnd->Rndm()*(GetNPixelsX()-1.) + 1.);
685 } while(ix<=0 || ix>GetNPixelsX());
687 iz=(Int_t)(rnd->Rndm()*(GetNPixelsZ()-1.) + 1.);
688 } while(iz<=0 || iz>GetNPixelsZ());
690 fMapA2->SetHit(iz,ix,signal);
691 } // end loop on masked pixels
693 //______________________________________________________________________
694 void AliITSsimulationSPD::CreateHistograms() {
698 fHis=new TObjArray(GetNPixelsZ());
699 for(i=0;i<GetNPixelsZ();i++) {
700 TString spdname("spd_");
702 sprintf(candnum,"%d",i+1);
703 spdname.Append(candnum);
704 (*fHis)[i] = new TH1F(spdname.Data(),"SPD maps",
705 GetNPixelsX(),0.,(Float_t) GetNPixelsX());
708 //______________________________________________________________________
709 void AliITSsimulationSPD::ResetHistograms() {
710 // Reset histograms for this detector
713 for(i=0;i<GetNPixelsZ();i++ ) {
714 if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
717 //______________________________________________________________________
718 void AliITSsimulationSPD::WriteSDigits(AliITSpList *pList){
719 // Fills the Summable digits Tree
721 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
723 pList->GetMaxMapIndex(ni,nj);
724 for(i=0;i<ni;i++)for(j=0;j<nj;j++){
725 if(pList->GetSignalOnly(i,j)>0.0){
726 aliITS->AddSumDigit(*(pList->GetpListItem(i,j)));
728 cout << "pListSPD: " << *(pList->GetpListItem(i,j)) << endl;
729 cout << " CreateSPDSDigit mod=" << fModule << " r,c=";
730 cout << i <<","<< j << " sig=" << fpList->GetSignalOnly(i,j);
731 cout << " noise=" << fpList->GetNoise(i,j) <<endl;
737 //______________________________________________________________________
738 void AliITSsimulationSPD::FillMapFrompList(AliITSpList *pList){
739 // Fills fMap2A from the pList of Summable digits
742 for(iz=0;iz<GetNPixelsZ();iz++)for(ix=0;ix<GetNPixelsX();ix++)
743 fMapA2->AddSignal(iz,ix,pList->GetSignal(iz,ix));