11 #include "AliITShit.h"
12 #include "AliITSdigit.h"
13 #include "AliITSmodule.h"
14 #include "AliITSMapA2.h"
15 #include "AliITSsimulationSPD.h"
16 #include "AliITSsegmentation.h"
17 #include "AliITSresponse.h"
22 ClassImp(AliITSsimulationSPD)
23 ////////////////////////////////////////////////////////////////////////
25 // Written by Boris Batyunya
28 // AliITSsimulationSPD is the simulation of SPDs
29 //________________________________________________________________________
32 AliITSsimulationSPD::AliITSsimulationSPD()
46 //_____________________________________________________________________________
48 AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) {
49 // standard constructor
55 fResponse->GetNoiseParam(fNoise,fBaseline);
57 fMapA2 = new AliITSMapA2(fSegmentation);
60 fNPixelsZ=fSegmentation->Npz();
61 fNPixelsX=fSegmentation->Npx();
65 //_____________________________________________________________________________
67 AliITSsimulationSPD::~AliITSsimulationSPD() {
79 //__________________________________________________________________________
80 AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
82 if(&source == this) return;
83 this->fMapA2 = source.fMapA2;
84 this->fNoise = source.fNoise;
85 this->fBaseline = source.fBaseline;
86 this->fNPixelsX = source.fNPixelsX;
87 this->fNPixelsZ = source.fNPixelsZ;
88 this->fHis = source.fHis;
92 //_________________________________________________________________________
94 AliITSsimulationSPD::operator=(const AliITSsimulationSPD &source) {
95 // Assignment operator
96 if(&source == this) return *this;
97 this->fMapA2 = source.fMapA2;
98 this->fNoise = source.fNoise;
99 this->fBaseline = source.fBaseline;
100 this->fNPixelsX = source.fNPixelsX;
101 this->fNPixelsZ = source.fNPixelsZ;
102 this->fHis = source.fHis;
105 //_____________________________________________________________________________
107 void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy)
111 const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
113 const Float_t kconv = 10000.; // cm -> microns
115 Float_t spdLength = fSegmentation->Dz();
116 Float_t spdWidth = fSegmentation->Dx();
118 Float_t difCoef, dum;
119 fResponse->DiffCoeff(difCoef,dum);
121 Float_t zPix0 = 1e+6;
122 Float_t xPix0 = 1e+6;
123 Float_t yPix0 = 1e+6;
124 Float_t yPrev = 1e+6;
128 Float_t zPitch = fSegmentation->Dpz(0);
129 Float_t xPitch = fSegmentation->Dpx(0);
131 //cout << "pitch per z: " << zPitch << endl;
132 //cout << "pitch per r*phi: " << xPitch << endl;
134 TObjArray *fHits = mod->GetHits();
135 Int_t nhits = fHits->GetEntriesFast();
137 // cout << "module, nhits ="<<module<<","<<nhits<< endl;
139 // Array of pointers to the label-signal list
141 Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;;
142 Float_t **pList = new Float_t* [maxNDigits];
143 memset(pList,0,sizeof(Float_t*)*maxNDigits);
144 Int_t indexRange[4] = {0,0,0,0};
146 // Fill detector maps with GEANT hits
147 // loop over hits in the module
149 Int_t lasttrack=-2,idhit=-1;
150 Int_t hit, iZi, jz, jx;
151 for (hit=0;hit<nhits;hit++) {
152 AliITShit *iHit = (AliITShit*) fHits->At(hit);
153 Int_t layer = iHit->GetLayer();
155 // work with the idtrack=entry number in the TreeH
156 //Int_t idhit,idtrack;
157 //mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);
158 //Int_t idtrack=mod->GetHitTrackIndex(hit);
159 // or store straight away the particle position in the array
162 if(iHit->StatusEntering()) idhit=hit;
163 Int_t itrack = iHit->GetTrack();
166 if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
168 // Int_t parent = iHit->GetParticle()->GetFirstMother();
169 Int_t partcode = iHit->GetParticle()->GetPdgCode();
171 // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
172 // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
174 Float_t px = iHit->GetPXL();
175 Float_t py = iHit->GetPYL();
176 Float_t pz = iHit->GetPZL();
177 Float_t pmod = 1000*sqrt(px*px+py*py+pz*pz);
180 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
184 // Get hit z and x(r*phi) cordinates for each module (detector)
187 Float_t zPix = kconv*iHit->GetZL();
188 Float_t xPix = kconv*iHit->GetXL();
189 Float_t yPix = kconv*iHit->GetYL();
192 Int_t status = iHit->GetTrackStatus();
195 if(TMath::Abs(zPix) > spdLength/2.) {
196 printf("!! Zpix outside = %f\n",zPix);
197 if(status == 66) zP0=100;
202 if (TMath::Abs(xPix) > spdWidth/2.) {
203 printf("!! Xpix outside = %f\n",xPix);
204 if (status == 66) xP0=100;
208 Float_t zP = (zPix + spdLength/2.)/1000.;
209 Float_t xP = (xPix + spdWidth/2.)/1000.;
213 // enter Si or after event in Si
220 if (layer == 1 && status == 66 && yPix > 71.) {
226 if (layer == 2 && status == 66 && yPix < -71.) {
231 Float_t depEnergy = iHit->GetIonization();
232 // skip if the input point to Si
233 if(depEnergy <= 0.) continue;
234 // skip if the input point is outside of Si, but the next
235 // point is inside of Si
236 if(zP0 > 90 || xP0 > 90) continue;
237 // if track returns to the opposite direction:
238 if (layer == 1 && yPix > yPrev) {
242 if (layer == 2 && yPix < yPrev) {
247 // take into account the holes diffusion inside the Silicon
248 // the straight line between the entrance and exit points in Si is
249 // divided into the several steps; the diffusion is considered
250 // for each end point of step and charge
251 // is distributed between the pixels through the diffusion.
254 // ---------- the diffusion in Z (beam) direction -------
256 Float_t charge = depEnergy*kEnToEl; // charge in e-
259 Float_t sigmaDif = 0.;
260 Float_t zdif = zPix - zPix0;
261 Float_t xdif = xPix - xPix0;
262 Float_t ydif = yPix - yPix0;
264 if(TMath::Abs(ydif) < 0.1) continue; // Ydif is not zero
266 Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
267 Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
268 Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
270 // number of the steps along the track:
272 if(ndX > ndZ) nsteps = ndX;
273 if(nsteps < 6) nsteps = 6; // minimum number of the steps
275 if(TMath::Abs(projDif) > 5.0) tang = ydif/projDif;
276 Float_t dCharge = charge/nsteps; // charge in e- for one step
277 Float_t dZ = zdif/nsteps;
278 Float_t dX = xdif/nsteps;
280 if (TMath::Abs(projDif) < 5.0 ) {
282 drPath = TMath::Abs(drPath); // drift path in cm
283 sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm
286 for (iZi = 1;iZi <= nsteps;iZi++) {
287 Float_t dZn = iZi*dZ;
288 Float_t dXn = iZi*dX;
289 Float_t zPixn = zPix0 + dZn;
290 Float_t xPixn = xPix0 + dXn;
292 if(TMath::Abs(projDif) >= 5.) {
293 Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
295 drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
296 drPath = TMath::Abs(drPath);
299 Float_t dProjn = projDif/nsteps;
300 drPath = (projDif-(iZi-1)*dProjn)*tang*1.e-4;
301 drPath = TMath::Abs(drPath);
303 sigmaDif = difCoef*sqrt(drPath);
304 sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
306 zPixn = (zPixn + spdLength/2.);
307 xPixn = (xPixn + spdWidth/2.);
309 fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
310 zPitch = fSegmentation->Dpz(nZpix);
311 fSegmentation->GetPadTxz(xPixn,zPixn);
312 // set the window for the integration
315 if(nZpix == 1) jzmin =2;
316 if(nZpix == fNPixelsZ) jzmax = 2;
320 if(nXpix == 1) jxmin =2;
321 if(nXpix == fNPixelsX) jxmax = 2;
323 Float_t zpix = nZpix;
324 Float_t dZright = zPitch*(zpix - zPixn);
325 Float_t dZleft = zPitch - dZright;
327 Float_t xpix = nXpix;
328 Float_t dXright = xPitch*(xpix - xPixn);
329 Float_t dXleft = xPitch - dXright;
336 for(jz=jzmin; jz <=jzmax; jz++) {
338 dZprev = -zPitch - dZleft;
347 dZnext = dZright + zPitch;
349 // kz changes from 1 to the fNofPixels(270)
350 Int_t kz = nZpix + jz -2;
352 Float_t zArg1 = dZprev/sigmaDif;
353 Float_t zArg2 = dZnext/sigmaDif;
354 Float_t zProb1 = TMath::Erfc(zArg1);
355 Float_t zProb2 = TMath::Erfc(zArg2);
356 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
359 // ----------- holes diffusion in X(r*phi) direction --------
362 for(jx=jxmin; jx <=jxmax; jx++) {
364 dXprev = -xPitch - dXleft;
373 dXnext = dXright + xPitch;
375 Int_t kx = nXpix + jx -2;
377 Float_t xArg1 = dXprev/sigmaDif;
378 Float_t xArg2 = dXnext/sigmaDif;
379 Float_t xProb1 = TMath::Erfc(xArg1);
380 Float_t xProb2 = TMath::Erfc(xArg2);
381 Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
387 indexRange[0]=indexRange[1]=index;
388 indexRange[2]=indexRange[3]=kx-1;
392 indexRange[0]=TMath::Min(indexRange[0],kz-1);
393 indexRange[1]=TMath::Max(indexRange[1],kz-1);
394 indexRange[2]=TMath::Min(indexRange[2],kx-1);
395 indexRange[3]=TMath::Max(indexRange[3],kx-1);
397 // build the list of digits for this module
398 Double_t signal=fMapA2->GetSignal(index,kx-1);
400 fMapA2->SetHit(index,kx-1,(double)signal);
407 if (status == 65) { // the step is inside of Si
414 GetList(itrack,idhit,pList,indexRange);
418 } // hit loop inside the module
421 // introduce the electronics effects and do zero-suppression
422 ChargeToSignal(pList);
431 //---------------------------------------------
432 void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
434 // lop over nonzero digits
438 for(int k=0;k<4;k++) {
439 if (indexRange[k] < 0) indexRange[k]=0;
442 for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
443 for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
445 Float_t signal=fMapA2->GetSignal(iz,ix);
447 if (!signal) continue;
449 Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
450 if(!pList[globalIndex]){
453 // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
456 pList[globalIndex] = new Float_t [9];
460 *pList[globalIndex] = -3.;
461 *(pList[globalIndex]+1) = -3.;
462 *(pList[globalIndex]+2) = -3.;
463 *(pList[globalIndex]+3) = 0.;
464 *(pList[globalIndex]+4) = 0.;
465 *(pList[globalIndex]+5) = 0.;
466 *(pList[globalIndex]+6) = -1.;
467 *(pList[globalIndex]+7) = -1.;
468 *(pList[globalIndex]+8) = -1.;
471 *pList[globalIndex] = (float)label;
472 *(pList[globalIndex]+3) = signal;
473 *(pList[globalIndex]+6) = (float)idhit;
477 // check the signal magnitude
479 Float_t highest = *(pList[globalIndex]+3);
480 Float_t middle = *(pList[globalIndex]+4);
481 Float_t lowest = *(pList[globalIndex]+5);
483 signal -= (highest+middle+lowest);
486 // compare the new signal with already existing list
489 if(signal<lowest) continue; // neglect this track
492 *(pList[globalIndex]+5) = middle;
493 *(pList[globalIndex]+4) = highest;
494 *(pList[globalIndex]+3) = signal;
496 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
497 *(pList[globalIndex]+1) = *pList[globalIndex];
498 *pList[globalIndex] = label;
500 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
501 *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
502 *(pList[globalIndex]+6) = idhit;
504 else if (signal>middle){
505 *(pList[globalIndex]+5) = middle;
506 *(pList[globalIndex]+4) = signal;
508 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
509 *(pList[globalIndex]+1) = label;
511 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
512 *(pList[globalIndex]+7) = idhit;
515 *(pList[globalIndex]+5) = signal;
516 *(pList[globalIndex]+2) = label;
517 *(pList[globalIndex]+8) = idhit;
520 } // end of loop pixels in x
521 } // end of loop over pixels in z
527 //---------------------------------------------
528 void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
530 // add noise and electronics, perform the zero suppression and add the
533 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
536 TRandom *random = new TRandom();
537 Float_t threshold = (float)fResponse->MinVal();
539 Int_t digits[3], tracks[3], hits[3],gi,j1;
543 for(Int_t iz=0;iz<fNPixelsZ;iz++){
544 for(Int_t ix=0;ix<fNPixelsX;ix++){
545 electronics = fBaseline + fNoise*random->Gaus();
546 signal = (float)fMapA2->GetSignal(iz,ix);
547 signal += electronics;
548 gi =iz*fNPixelsX+ix; // global index
549 if (signal > threshold) {
556 tracks[j1] = (Int_t)(*(pList[gi]+j1));
557 hits[j1] = (Int_t)(*(pList[gi]+j1+6));
559 tracks[j1]=-2; //noise
565 if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
571 if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
575 if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
579 if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
585 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
587 if(pList[gi]) delete [] pList[gi];
595 //____________________________________________
597 void AliITSsimulationSPD::CreateHistograms()
599 // create 1D histograms for tests
601 printf("SPD - create histograms\n");
603 fHis=new TObjArray(fNPixelsZ);
604 TString spdName("spd_");
605 for (Int_t i=0;i<fNPixelsZ;i++) {
607 sprintf(pixelz,"%d",i+1);
608 spdName.Append(pixelz);
609 (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
610 fNPixelsX,0.,(Float_t) fNPixelsX);
614 //____________________________________________
616 void AliITSsimulationSPD::ResetHistograms()
619 // Reset histograms for this detector
621 for ( int i=0;i<fNPixelsZ;i++ ) {
622 if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();