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 <Riostream.h>
23 #include <TParticle.h>
27 #include "AliITShit.h"
28 #include "AliITSdigitSPD.h"
29 #include "AliITSmodule.h"
30 #include "AliITSMapA2.h"
31 #include "AliITSpList.h"
32 #include "AliITSsimulationSPD.h"
33 #include "AliITSsegmentation.h"
34 #include "AliITSresponse.h"
35 #include "AliITSsegmentationSPD.h"
36 #include "AliITSresponseSPD.h"
40 ClassImp(AliITSsimulationSPD)
41 ////////////////////////////////////////////////////////////////////////
43 // Written by Rocco Caliandro
44 // from a model developed with T. Virgili and R.A. Fini
47 // AliITSsimulationSPD is the simulation of SPDs
49 //______________________________________________________________________
50 AliITSsimulationSPD::AliITSsimulationSPD(){
51 // Default constructor
64 //______________________________________________________________________
65 AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg,
66 AliITSresponse *resp) {
67 // Standard constructor
80 Init((AliITSsegmentationSPD*)seg,(AliITSresponseSPD*)resp);
82 //______________________________________________________________________
83 void AliITSsimulationSPD::Init(AliITSsegmentationSPD *seg,
84 AliITSresponseSPD *resp) {
85 // Initilizes the variables of AliITSsimulation SPD.
90 fMapA2 = new AliITSMapA2(fSegmentation);
91 fpList = new AliITSpList(GetNPixelsZ()+1,GetNPixelsX()+1);
93 fResponse->Thresholds(fThresh,fSigma);
94 fResponse->GetNoiseParam(fCouplCol,fCouplRow);
95 fNPixelsZ = fSegmentation->Npz();
96 fNPixelsX = fSegmentation->Npx();
99 //______________________________________________________________________
100 AliITSsimulationSPD::~AliITSsimulationSPD() {
111 //______________________________________________________________________
112 AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source) :
113 AliITSsimulation(source){
116 if(&source == this) return;
118 this->fMapA2 = source.fMapA2;
119 this->fHis = source.fHis;
121 this->fThresh = source.fThresh;
122 this->fSigma = source.fSigma;
123 this->fCouplCol = source.fCouplCol;
124 this->fCouplRow = source.fCouplRow;
125 this->fNPixelsX = source.fNPixelsX;
126 this->fNPixelsZ = source.fNPixelsZ;
130 //______________________________________________________________________
131 AliITSsimulationSPD& AliITSsimulationSPD::operator=(const AliITSsimulationSPD
133 // Assignment operator
135 if(&source == this) return *this;
137 this->fMapA2 = source.fMapA2;
138 this->fHis = source.fHis;
140 this->fThresh = source.fThresh;
141 this->fSigma = source.fSigma;
142 this->fCouplCol = source.fCouplCol;
143 this->fCouplRow = source.fCouplRow;
144 this->fNPixelsX = source.fNPixelsX;
145 this->fNPixelsZ = source.fNPixelsZ;
149 //______________________________________________________________________
150 void AliITSsimulationSPD::InitSimulationModule(Int_t module,Int_t event){
151 // Creates maps to build the list of tracks for each sumable digit
153 // Int_t module // Module number to be simulated
154 // Int_t event // Event number to be simulated
165 //______________________________________________________________________
166 void AliITSsimulationSPD::FinishSDigitiseModule(){
167 // Does the Sdigits to Digits work
175 SDigitsToDigits(fModule,fpList);
177 //______________________________________________________________________
178 void AliITSsimulationSPD::SDigitiseModule(AliITSmodule *mod, Int_t dummy0,
180 // Sum digitize module
181 if (!(mod->GetNhits())) return; // if module has no hits then no Sdigits.
182 Int_t number = 10000;
183 Int_t *frowpixel = new Int_t[number];
184 Int_t *fcolpixel = new Int_t[number];
185 Double_t *fenepixel = new Double_t[number];
187 dummy0 = dummy1; // remove unsued variable warning.
188 fModule = mod->GetIndex();
190 // Array of pointers to store the track index of the digits
191 // leave +1, otherwise pList crashes when col=256, row=192
193 HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,fpList);
195 WriteSDigits(fpList);
204 //______________________________________________________________________
205 void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t dummy0,
207 // digitize module. Also need to digitize modules with only noise.
209 Int_t number = 10000;
210 Int_t *frowpixel = new Int_t[number];
211 Int_t *fcolpixel = new Int_t[number];
212 Double_t *fenepixel = new Double_t[number];
214 dummy0 = dummy1; // remove unsued variable warning.
215 // Array of pointers to store the track index of the digits
216 // leave +1, otherwise pList crashes when col=256, row=192
217 fModule = mod->GetIndex();
219 SetFluctuations(fpList,fModule);
221 HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,fpList);
223 // apply mask to SPD module
226 CreateDigit(fModule,fpList);
235 //______________________________________________________________________
236 void AliITSsimulationSPD::SDigitsToDigits(Int_t module,AliITSpList *pList) {
237 // sum digits to Digits.
240 cout << "Entering AliITSsimulatinSPD::SDigitsToDigits for module=";
241 cout << module << endl;
246 SetFluctuations(pList,module);
248 fMapA2->ClearMap(); // since noise is in pList aready. Zero Map so that
249 // noise is not doubled when calling FillMapFrompList.
251 FillMapFrompList(pList);
253 // apply mask to SPD module
256 CreateDigit(module,pList);
261 //______________________________________________________________________
262 void AliITSsimulationSPD::UpdateMapSignal(Int_t row,Int_t col,Int_t trk,
263 Int_t hit,Int_t mod,Double_t ene,
264 AliITSpList *pList) {
265 // updates the Map of signal, adding the energy (ene) released by
268 fMapA2->AddSignal(row,col,ene);
269 pList->AddSignal(row,col,trk,hit,mod,ene);
271 //______________________________________________________________________
272 void AliITSsimulationSPD::UpdateMapNoise(Int_t row,Int_t col,Int_t mod,
273 Double_t ene,AliITSpList *pList) {
274 // updates the Map of noise, adding the energy (ene) give my noise
276 fMapA2->AddSignal(row,col,ene);
277 pList->AddNoise(row,col,mod,ene);
279 //______________________________________________________________________
280 void AliITSsimulationSPD::HitsToAnalogDigits(AliITSmodule *mod,
281 Int_t *frowpixel,Int_t *fcolpixel,
283 AliITSpList *pList) {
284 // Loops over all hits to produce Analog/floting point digits. This
285 // is also the first task in producing standard digits.
287 // loop over hits in the module
288 Int_t hitpos,nhits = mod->GetNhits();
289 for (hitpos=0;hitpos<nhits;hitpos++) {
290 HitToDigit(mod,hitpos,frowpixel,fcolpixel,fenepixel,pList);
291 }// end loop over digits
293 //______________________________________________________________________
294 void AliITSsimulationSPD::HitToDigit(AliITSmodule *mod,Int_t hitpos,
295 Int_t *frowpixel,Int_t *fcolpixel,
296 Double_t *fenepixel,AliITSpList *pList) {
297 // Steering function to determine the digits associated to a given
299 // The digits are created by charge sharing (ChargeSharing) and by
300 // capacitive coupling (SetCoupling). At all the created digits is
301 // associated the track number of the hit (ntrack)
302 Double_t x1l=0.0,y1l=0.0,z1l=0.0,x2l=0.0,y2l=0.0,z2l=0.0;
303 Int_t r1,r2,c1,c2,row,col,npixel = 0;
305 Double_t ene=0.0,etot=0.0;
306 const Float_t kconv = 10000.; // cm -> microns
307 const Float_t kconv1= 0.277e9; // GeV -> electrons equivalent
309 if(!(mod->LineSegmentL(hitpos,x1l,x2l,y1l,y2l,z1l,z2l,etot,ntrack)))return;
311 x2l += x1l; y2l += y1l; z2l += z1l; // Convert to ending coordinate.
312 // positions shifted and converted in microns
313 x1l = x1l*kconv + fSegmentation->Dx()/2.;
314 z1l = z1l*kconv + fSegmentation->Dz()/2.;
315 // positions shifted and converted in microns
316 x2l = x2l*kconv + fSegmentation->Dx()/2.;
317 z2l = z2l*kconv + fSegmentation->Dz()/2.;
318 etot *= kconv1; // convert from GeV to electrons equivalent.
319 Int_t module = mod->GetIndex();
321 // to account for the effective sensitive area
322 // introduced in geometry
323 if (z1l<0 || z1l>fSegmentation->Dz()) return;
324 if (z2l<0 || z2l>fSegmentation->Dz()) return;
325 if (x1l<0 || x1l>fSegmentation->Dx()) return;
326 if (x2l<0 || x2l>fSegmentation->Dx()) return;
328 //Get the col and row number starting from 1
329 // the x direction is not inverted for the second layer!!!
330 fSegmentation->GetPadIxz(x1l, z1l, c1, r1);
331 fSegmentation->GetPadIxz(x2l, z2l, c2, r2);
333 // to account for unexpected equal entrance and
335 if (x1l==x2l) x2l=x2l+x2l*0.1;
336 if (z1l==z2l) z2l=z2l+z2l*0.1;
338 if ((r1==r2) && (c1==c2)){
341 frowpixel[npixel-1] = r1;
342 fcolpixel[npixel-1] = c1;
343 fenepixel[npixel-1] = etot;
346 ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
347 npixel,frowpixel,fcolpixel,fenepixel);
348 } // end if r1==r2 && c1==c2.
350 for (Int_t npix=0;npix<npixel;npix++){
351 row = frowpixel[npix];
352 col = fcolpixel[npix];
353 ene = fenepixel[npix];
354 UpdateMapSignal(row,col,ntrack,hitpos,module,ene,pList);
355 // Starting capacitive coupling effect
356 SetCoupling(row,col,ntrack,hitpos,module,pList);
359 //______________________________________________________________________
360 void AliITSsimulationSPD::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
361 Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
362 Int_t r2,Float_t etot,
363 Int_t &npixel,Int_t *frowpixel,
364 Int_t *fcolpixel,Double_t *fenepixel){
365 // Take into account the geometrical charge sharing when the track
366 // crosses more than one pixel.
370 <img src="picts/ITS/barimodel_2.gif">
373 <font size=+2 color=red>
374 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
380 Float_t xa,za,xb,zb,dx,dz,dtot,refr,refm,refc;
382 Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
383 Int_t dirx,dirz,rb,cb;
384 Int_t flag,flagrow,flagcol;
394 dtot = TMath::Sqrt((dx*dx)+(dz*dz));
395 if (dtot==0.0) dtot = 0.01;
396 dirx = (Int_t) TMath::Sign((Float_t)1,dx);
397 dirz = (Int_t) TMath::Sign((Float_t)1,dz);
399 // calculate the x coordinate of the pixel in the next column
400 // and the z coordinate of the pixel in the next row
403 fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos);
405 Float_t xsize = fSegmentation->Dpx(0);
406 Float_t zsize = fSegmentation->Dpz(r1-1);
408 if (dirx == 1) refr = xpos+xsize/2.;
409 else refr = xpos-xsize/2.;
411 if (dirz == 1) refn = zpos+zsize/2.;
412 else refn = zpos-zsize/2.;
418 // calculate the x coordinate of the intersection with the pixel
419 // in the next cell in row direction
421 refm = dx*((refn - z1l)/dz) + x1l;
423 refm = refr+dirx*xsize;
425 // calculate the z coordinate of the intersection with the pixel
426 // in the next cell in column direction
428 refc = dz*((refr - x1l)/dx) + z1l;
430 refc = refn+dirz*zsize;
437 if ((arefm < arefr) && (arefn < arefc)){
438 // the track goes in the pixel in the next cell in row direction
445 if (rb == r2) flagrow=1;
450 // shift to the pixel in the next cell in row direction
451 Float_t zsizeNext = fSegmentation->Dpz(rb-1);
452 //to account for cell at the borders of the detector
453 if(zsizeNext==0) zsizeNext = zsize;
454 refn += zsizeNext*dirz;
456 // the track goes in the pixel in the next cell in column direction
463 if (cb == c2) flagcol=1;
467 } // end ifaxb > ax2l
469 // shift to the pixel in the next cell in column direction
470 Float_t xsizeNext = fSegmentation->Dpx(cb-1);
471 //to account for cell at the borders of the detector
472 if(xsizeNext==0) xsizeNext = xsize;
473 refr += xsizeNext*dirx;
474 } // end if (arefm < arefr) && (arefn < arefc)
476 //calculate the energy lost in the crossed pixel
477 epar = TMath::Sqrt((xb-xa)*(xb-xa)+(zb-za)*(zb-za));
478 epar = etot*(epar/dtot);
480 //store row, column and energy lost in the crossed pixel
481 frowpixel[npixel] = r1;
482 fcolpixel[npixel] = c1;
483 fenepixel[npixel] = epar;
486 // the exit point of the track is reached
487 if (epar == 0) flag = 1;
488 if ((r1 == r2) && (c1 == c2)) flag = 1;
497 //______________________________________________________________________
498 void AliITSsimulationSPD::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
499 Int_t idhit,Int_t module,
500 AliITSpList *pList) {
501 // Take into account the coupling between adiacent pixels.
502 // The parameters probcol and probrow are the probability of the
503 // signal in one pixel shared in the two adjacent pixels along
504 // the column and row direction, respectively.
508 <img src="picts/ITS/barimodel_3.gif">
511 <font size=+2 color=red>
512 <a href="mailto:tiziano.virgili@cern.ch"></a>.
518 Double_t pulse1,pulse2;
519 Float_t couplR=0.0,couplC=0.0;
522 GetCouplings(couplR,couplC);
525 pulse1 = fMapA2->GetSignal(row,col);
527 for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
531 xr = gRandom->Rndm();
532 // if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
533 if ((j1<0) || (j1>GetNPixelsZ()-1) || (xr>couplR)){
537 UpdateMapSignal(j1,col,ntrack,idhit,module,pulse1,pList);
539 flag = 1; // only first next!!
542 // loop in column direction
546 xr = gRandom->Rndm();
547 // if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
548 if ((j2<0) || (j2>GetNPixelsX()-1) || (xr>couplC)){
552 UpdateMapSignal(row,j2,ntrack,idhit,module,pulse2,pList);
554 flag = 1; // only first next!!
559 //______________________________________________________________________
560 void AliITSsimulationSPD::SetCouplingOld(Int_t row, Int_t col, Int_t ntrack,
561 Int_t idhit,Int_t module,
562 AliITSpList *pList) {
563 // Take into account the coupling between adiacent pixels.
564 // The parameters probcol and probrow are the fractions of the
565 // signal in one pixel shared in the two adjacent pixels along
566 // the column and row direction, respectively.
570 <img src="picts/ITS/barimodel_3.gif">
573 <font size=+2 color=red>
574 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
580 Double_t pulse1,pulse2;
581 Float_t couplR=0.0,couplC=0.0;
583 GetCouplings(couplR,couplC);
586 pulse1 = fMapA2->GetSignal(row,col);
588 for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
592 if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
593 pulse1 = fMapA2->GetSignal(row,col);
597 UpdateMapSignal(j1,col,ntrack,idhit,module,pulse1,pList);
601 // loop in column direction
605 if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
606 pulse2 = fMapA2->GetSignal(row,col);
610 UpdateMapSignal(row,j2,ntrack,idhit,module,pulse2,pList);
616 //______________________________________________________________________
617 void AliITSsimulationSPD::CreateDigit(Int_t module,AliITSpList *pList) {
618 // The pixels are fired if the energy deposited inside them is above
619 // the threshold parameter ethr. Fired pixed are interpreted as digits
620 // and stored in the file digitfilename. One also needs to write out
621 // cases when there is only noise (nhits==0).
623 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
625 Int_t size = AliITSdigitSPD::GetNTracks();
626 Int_t * digits = new Int_t[size];
627 Int_t * tracks = new Int_t[size];
628 Int_t * hits = new Int_t[size];
629 Float_t * charges = new Float_t[size];
632 module=0; // remove unused variable warning.
633 for(j1=0;j1<size;j1++){tracks[j1]=-3;hits[j1]=-1;charges[j1]=0.0;}
634 for (Int_t r=1;r<=GetNPixelsZ();r++) {
635 for (Int_t c=1;c<=GetNPixelsX();c++) {
636 // check if the deposited energy in a pixel is above the
638 Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
639 if ( signal > GetThreshold()) {
640 digits[0] = r-1; // digits starts from 0
641 digits[1] = c-1; // digits starts from 0
643 digits[2] = (Int_t) signal; // the signal is stored in
645 for(j1=0;j1<size;j1++){
646 if(j1<pList->GetNEnteries()){
647 tracks[j1] = pList->GetTrack(r,c,j1);
648 hits[j1] = pList->GetHit(r,c,j1);
656 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
658 cout << " CreateSPDDigit mod=" << fModule << " r,c=" << r;
659 cout <<","<<c<< " sig=" << fpList->GetSignalOnly(r,c);
660 cout << " noise=" << fpList->GetNoise(r,c);
661 cout << " Msig="<< signal << " Thres=" << GetThreshold()<<endl;
663 } // end if of threshold condition
672 //______________________________________________________________________
673 void AliITSsimulationSPD::SetFluctuations(AliITSpList *pList,Int_t module) {
674 // Set the electronic noise and threshold non-uniformities to all the
675 // pixels in a detector.
676 // The parameter fSigma is the squared sum of the sigma due to noise
677 // and the sigma of the threshold distribution among pixels.
681 <img src="picts/ITS/barimodel_1.gif">
684 <font size=+2 color=red>
685 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
690 Float_t thr=0.0,sigm=0.0;
691 Double_t signal,sigma;
694 GetThresholds(thr,sigm);
695 sigma = (Double_t) sigm;
696 for(iz=1;iz<=GetNPixelsZ();iz++){
697 for(ix=1;ix<=GetNPixelsX();ix++){
698 signal = sigma*gRandom->Gaus();
699 fMapA2->SetHit(iz,ix,signal);
700 // insert in the label-signal-hit list the pixels fired
702 pList->AddNoise(iz,ix,module,signal);
703 } // end of loop on pixels
704 } // end of loop on pixels
706 //______________________________________________________________________
707 void AliITSsimulationSPD::SetMask() {
708 // Apply a mask to the SPD module. 1% of the pixel channels are
709 // masked. When the database will be ready, the masked pixels
710 // should be read from it.
714 Float_t perc = ((AliITSresponseSPD*)fResponse)->GetFractionDead();
715 // in this way we get the same set of random numbers for all runs.
716 // This is a cluge for now.
717 static TRandom *rnd = new TRandom();
719 totMask= perc*GetNPixelsZ()*GetNPixelsX();
720 for(im=1;im<totMask;im++){
722 ix=(Int_t)(rnd->Rndm()*(GetNPixelsX()-1.) + 1.);
723 } while(ix<=0 || ix>GetNPixelsX());
725 iz=(Int_t)(rnd->Rndm()*(GetNPixelsZ()-1.) + 1.);
726 } while(iz<=0 || iz>GetNPixelsZ());
728 fMapA2->SetHit(iz,ix,signal);
729 } // end loop on masked pixels
731 //______________________________________________________________________
732 void AliITSsimulationSPD::CreateHistograms() {
736 fHis=new TObjArray(GetNPixelsZ());
737 for(i=0;i<GetNPixelsZ();i++) {
738 TString spdname("spd_");
740 sprintf(candnum,"%d",i+1);
741 spdname.Append(candnum);
742 (*fHis)[i] = new TH1F(spdname.Data(),"SPD maps",
743 GetNPixelsX(),0.,(Float_t) GetNPixelsX());
746 //______________________________________________________________________
747 void AliITSsimulationSPD::ResetHistograms() {
748 // Reset histograms for this detector
751 for(i=0;i<GetNPixelsZ();i++ ) {
752 if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
755 //______________________________________________________________________
756 void AliITSsimulationSPD::WriteSDigits(AliITSpList *pList){
757 // Fills the Summable digits Tree
759 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
761 pList->GetMaxMapIndex(ni,nj);
762 for(i=0;i<ni;i++)for(j=0;j<nj;j++){
763 if(pList->GetSignalOnly(i,j)>0.0){
764 aliITS->AddSumDigit(*(pList->GetpListItem(i,j)));
766 cout << "pListSPD: " << *(pList->GetpListItem(i,j)) << endl;
767 cout << " CreateSPDSDigit mod=" << fModule << " r,c=";
768 cout << i <<","<< j << " sig=" << fpList->GetSignalOnly(i,j);
769 cout << " noise=" << fpList->GetNoise(i,j) <<endl;
775 //______________________________________________________________________
776 void AliITSsimulationSPD::FillMapFrompList(AliITSpList *pList){
777 // Fills fMap2A from the pList of Summable digits
780 for(iz=0;iz<GetNPixelsZ();iz++)for(ix=0;ix<GetNPixelsX();ix++)
781 fMapA2->AddSignal(iz,ix,pList->GetSignal(iz,ix));