11 #include "AliITShit.h"
12 #include "AliITSdigit.h"
13 #include "AliITSmodule.h"
14 #include "AliITSMapA2.h"
15 #include "AliITSsimulationSPDdubna.h"
16 #include "AliITSsegmentation.h"
17 #include "AliITSresponse.h"
22 ClassImp(AliITSsimulationSPDdubna)
23 ////////////////////////////////////////////////////////////////////////
25 // Written by Boris Batyunya
28 // AliITSsimulationSPDdubna is the simulation of SPDs
29 //________________________________________________________________________
32 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna()
46 //_____________________________________________________________________________
48 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(AliITSsegmentation *seg, AliITSresponse *resp) {
49 // standard constructor
55 fResponse->GetNoiseParam(fNoise,fBaseline);
57 fMapA2 = new AliITSMapA2(fSegmentation);
61 fNPixelsZ=fSegmentation->Npz();
62 fNPixelsX=fSegmentation->Npx();
66 //_____________________________________________________________________________
68 AliITSsimulationSPDdubna::~AliITSsimulationSPDdubna() {
80 //__________________________________________________________________________
81 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(const AliITSsimulationSPDdubna &source){
83 if(&source == this) return;
84 this->fMapA2 = source.fMapA2;
85 this->fNoise = source.fNoise;
86 this->fBaseline = source.fBaseline;
87 this->fNPixelsX = source.fNPixelsX;
88 this->fNPixelsZ = source.fNPixelsZ;
89 this->fHis = source.fHis;
93 //_________________________________________________________________________
94 AliITSsimulationSPDdubna&
95 AliITSsimulationSPDdubna::operator=(const AliITSsimulationSPDdubna &source) {
96 // Assignment operator
97 if(&source == this) return *this;
98 this->fMapA2 = source.fMapA2;
99 this->fNoise = source.fNoise;
100 this->fBaseline = source.fBaseline;
101 this->fNPixelsX = source.fNPixelsX;
102 this->fNPixelsZ = source.fNPixelsZ;
103 this->fHis = source.fHis;
106 //_____________________________________________________________________________
108 void AliITSsimulationSPDdubna::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy)
112 const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
114 const Float_t kconv = 10000.; // cm -> microns
116 Float_t spdLength = fSegmentation->Dz();
117 Float_t spdWidth = fSegmentation->Dx();
118 Float_t spdThickness = fSegmentation->Dy();
119 Float_t difCoef, dum;
120 fResponse->DiffCoeff(difCoef,dum);
121 if(spdThickness > 290) difCoef = 0.00613;
123 Float_t zPix0 = 1e+6;
124 Float_t xPix0 = 1e+6;
125 Float_t yPrev = 1e+6;
127 Float_t zPitch = fSegmentation->Dpz(0);
128 Float_t xPitch = fSegmentation->Dpx(0);
130 TObjArray *fHits = mod->GetHits();
131 Int_t nhits = fHits->GetEntriesFast();
134 cout<<"len,wid,thickness,nx,nz,pitchx,pitchz,difcoef ="<<spdLength<<","<<spdWidth<<","<<spdThickness<<","<<fNPixelsX<<","<<fNPixelsZ<<","<<xPitch<<","<<zPitch<<","<<difCoef<<endl;
135 // Array of pointers to the label-signal list
137 Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;;
138 Float_t **pList = new Float_t* [maxNDigits];
139 memset(pList,0,sizeof(Float_t*)*maxNDigits);
140 Int_t indexRange[4] = {0,0,0,0};
142 // Fill detector maps with GEANT hits
143 // loop over hits in the module
146 Int_t hit, iZi, jz, jx;
148 cout<<"SPDdubna: module,nhits ="<<module<<","<<nhits<<endl;
149 for (hit=0;hit<nhits;hit++) {
150 AliITShit *iHit = (AliITShit*) fHits->At(hit);
151 //Int_t layer = iHit->GetLayer();
152 Float_t yPix0 = -spdThickness/2;
154 // work with the idtrack=entry number in the TreeH
155 //Int_t idhit,idtrack; //!
156 //mod->GetHitTrackAndHitIndex(hit,idtrack,idhit); //!
157 //Int_t idtrack=mod->GetHitTrackIndex(hit);
158 // or store straight away the particle position in the array
160 if(iHit->StatusEntering()) idhit=hit;
161 Int_t itrack = iHit->GetTrack();
164 if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
166 //Int_t parent = iHit->GetParticle()->GetFirstMother();
167 Int_t partcode = iHit->GetParticle()->GetPdgCode();
169 // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
170 // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
173 Float_t px = iHit->GetPXL(); // the momenta at the
174 Float_t py = iHit->GetPYL(); // each GEANT step
175 Float_t pz = iHit->GetPZL();
176 Float_t ptot = 1000*sqrt(px*px+py*py+pz*pz);
179 Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
184 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
188 // Get hit z and x(r*phi) cordinates for each module (detector)
191 Float_t zPix = kconv*iHit->GetZL();
192 Float_t xPix = kconv*iHit->GetXL();
193 Float_t yPix = kconv*iHit->GetYL();
196 Int_t status = iHit->GetTrackStatus();
199 if(zPix > spdLength/2) {
200 //cout<<"!!! SPD: z outside ="<<zPix<<endl;
201 zPix = spdLength/2 - 10;
203 if(zPix < 0 && zPix < -spdLength/2) {
204 //cout<<"!!! SPD: z outside ="<<zPix<<endl;
205 zPix = -spdLength/2 + 10;
207 if(xPix > spdWidth/2) {
208 //cout<<"!!! SPD: x outside ="<<xPix<<endl;
209 xPix = spdWidth/2 - 10;
211 if(xPix < 0 && xPix < -spdWidth/2) {
212 //cout<<"!!! SPD: x outside ="<<xPix<<endl;
213 xPix = -spdWidth/2 + 10;
217 // enter Si or after event in Si
224 Float_t depEnergy = iHit->GetIonization();
225 // skip if the input point to Si
227 if(depEnergy <= 0.) continue;
229 // if track returns to the opposite direction:
235 // take into account the holes diffusion inside the Silicon
236 // the straight line between the entrance and exit points in Si is
237 // divided into the several steps; the diffusion is considered
238 // for each end point of step and charge
239 // is distributed between the pixels through the diffusion.
242 // ---------- the diffusion in Z (beam) direction -------
244 Float_t charge = depEnergy*kEnToEl; // charge in e-
247 Float_t sigmaDif = 0.;
248 Float_t zdif = zPix - zPix0;
249 Float_t xdif = xPix - xPix0;
250 Float_t ydif = TMath::Abs(yPix - yPrev);
251 Float_t ydif0 = TMath::Abs(yPrev - yPix0);
253 if(ydif < 1) continue; // ydif is not zero
255 Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
257 Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
258 Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
260 // number of the steps along the track:
262 if(ndX > ndZ) nsteps = ndX;
263 if(nsteps < 20) nsteps = 20; // minimum number of the steps
266 drPath = (yPix-yPix0)*1.e-4;
267 drPath = TMath::Abs(drPath); // drift path in cm
268 sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm
269 sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
273 if(projDif > 5) tang = ydif/projDif;
274 Float_t dCharge = charge/nsteps; // charge in e- for one step
275 Float_t dZ = zdif/nsteps;
276 Float_t dX = xdif/nsteps;
278 for (iZi = 1;iZi <= nsteps;iZi++) {
279 Float_t dZn = iZi*dZ;
280 Float_t dXn = iZi*dX;
281 Float_t zPixn = zPix0 + dZn;
282 Float_t xPixn = xPix0 + dXn;
285 Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
286 drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
288 drPath = TMath::Abs(drPath) + ydif0*1.e-4;
291 drPath = ydif0*1.e-4 - TMath::Abs(drPath);
292 drPath = TMath::Abs(drPath);
294 sigmaDif = difCoef*sqrt(drPath);
295 sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
298 zPixn = (zPixn + spdLength/2.);
299 xPixn = (xPixn + spdWidth/2.);
301 fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
302 zPitch = fSegmentation->Dpz(nZpix);
303 fSegmentation->GetPadTxz(xPixn,zPixn);
304 // set the window for the integration
307 if(nZpix == 1) jzmin =2;
308 if(nZpix == fNPixelsZ) jzmax = 2;
312 if(nXpix == 1) jxmin =2;
313 if(nXpix == fNPixelsX) jxmax = 2;
315 Float_t zpix = nZpix;
316 Float_t dZright = zPitch*(zpix - zPixn);
317 Float_t dZleft = zPitch - dZright;
319 Float_t xpix = nXpix;
320 Float_t dXright = xPitch*(xpix - xPixn);
321 Float_t dXleft = xPitch - dXright;
328 for(jz=jzmin; jz <=jzmax; jz++) {
330 dZprev = -zPitch - dZleft;
339 dZnext = dZright + zPitch;
341 // kz changes from 1 to the fNofPixels(270)
342 Int_t kz = nZpix + jz -2;
344 Float_t zArg1 = dZprev/sigmaDif;
345 Float_t zArg2 = dZnext/sigmaDif;
346 Float_t zProb1 = TMath::Erfc(zArg1);
347 Float_t zProb2 = TMath::Erfc(zArg2);
348 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
351 // ----------- holes diffusion in X(r*phi) direction --------
354 for(jx=jxmin; jx <=jxmax; jx++) {
356 dXprev = -xPitch - dXleft;
365 dXnext = dXright + xPitch;
367 Int_t kx = nXpix + jx -2;
369 Float_t xArg1 = dXprev/sigmaDif;
370 Float_t xArg2 = dXnext/sigmaDif;
371 Float_t xProb1 = TMath::Erfc(xArg1);
372 Float_t xProb2 = TMath::Erfc(xArg2);
373 Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
379 indexRange[0]=indexRange[1]=index;
380 indexRange[2]=indexRange[3]=kx-1;
384 indexRange[0]=TMath::Min(indexRange[0],kz-1);
385 indexRange[1]=TMath::Max(indexRange[1],kz-1);
386 indexRange[2]=TMath::Min(indexRange[2],kx-1);
387 indexRange[3]=TMath::Max(indexRange[3],kx-1);
389 // build the list of digits for this module
390 Double_t signal=fMapA2->GetSignal(index,kx-1);
392 fMapA2->SetHit(index,kx-1,(double)signal);
399 if (status == 65) { // the step is inside of Si
406 GetList(itrack,idhit,pList,indexRange);
410 } // hit loop inside the module
413 // introduce the electronics effects and do zero-suppression
414 ChargeToSignal(pList);
423 //---------------------------------------------
424 void AliITSsimulationSPDdubna::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
426 // lop over nonzero digits
430 for(int k=0;k<4;k++) {
431 if (indexRange[k] < 0) indexRange[k]=0;
434 for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
435 for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
437 Float_t signal=fMapA2->GetSignal(iz,ix);
439 if (!signal) continue;
441 Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
442 if(!pList[globalIndex]){
445 // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
448 pList[globalIndex] = new Float_t [9];
452 *pList[globalIndex] = -3.;
453 *(pList[globalIndex]+1) = -3.;
454 *(pList[globalIndex]+2) = -3.;
455 *(pList[globalIndex]+3) = 0.;
456 *(pList[globalIndex]+4) = 0.;
457 *(pList[globalIndex]+5) = 0.;
458 *(pList[globalIndex]+6) = -1.;
459 *(pList[globalIndex]+7) = -1.;
460 *(pList[globalIndex]+8) = -1.;
463 *pList[globalIndex] = (float)label;
464 *(pList[globalIndex]+3) = signal;
465 *(pList[globalIndex]+6) = (float)idhit;
469 // check the signal magnitude
471 Float_t highest = *(pList[globalIndex]+3);
472 Float_t middle = *(pList[globalIndex]+4);
473 Float_t lowest = *(pList[globalIndex]+5);
475 signal -= (highest+middle+lowest);
478 // compare the new signal with already existing list
481 if(signal<lowest) continue; // neglect this track
484 *(pList[globalIndex]+5) = middle;
485 *(pList[globalIndex]+4) = highest;
486 *(pList[globalIndex]+3) = signal;
488 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
489 *(pList[globalIndex]+1) = *pList[globalIndex];
490 *pList[globalIndex] = label;
492 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
493 *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
494 *(pList[globalIndex]+6) = idhit;
496 else if (signal>middle){
497 *(pList[globalIndex]+5) = middle;
498 *(pList[globalIndex]+4) = signal;
500 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
501 *(pList[globalIndex]+1) = label;
503 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
504 *(pList[globalIndex]+7) = idhit;
507 *(pList[globalIndex]+5) = signal;
508 *(pList[globalIndex]+2) = label;
509 *(pList[globalIndex]+8) = idhit;
512 } // end of loop pixels in x
513 } // end of loop over pixels in z
519 //---------------------------------------------
520 void AliITSsimulationSPDdubna::ChargeToSignal(Float_t **pList)
522 // add noise and electronics, perform the zero suppression and add the
525 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
528 Float_t threshold = (float)fResponse->MinVal();
530 Int_t digits[3], tracks[3], hits[3],gi,j1;
534 for(Int_t iz=0;iz<fNPixelsZ;iz++){
535 for(Int_t ix=0;ix<fNPixelsX;ix++){
536 electronics = fBaseline + fNoise*gRandom->Gaus();
537 signal = (float)fMapA2->GetSignal(iz,ix);
538 signal += electronics;
539 gi =iz*fNPixelsX+ix; // global index
540 if (signal > threshold) {
546 //b.b. tracks[j1]=-3;
547 tracks[j1] = (Int_t)(*(pList[gi]+j1));
548 hits[j1] = (Int_t)(*(pList[gi]+j1+6));
550 tracks[j1]=-2; //noise
556 if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
562 if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
566 if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
570 if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
576 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
578 if(pList[gi]) delete [] pList[gi];
586 //____________________________________________
588 void AliITSsimulationSPDdubna::CreateHistograms()
590 // create 1D histograms for tests
592 printf("SPD - create histograms\n");
594 fHis=new TObjArray(fNPixelsZ);
595 TString spdName("spd_");
596 for (Int_t i=0;i<fNPixelsZ;i++) {
598 sprintf(pixelz,"%d",i+1);
599 spdName.Append(pixelz);
600 //PH (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
601 //PH fNPixelsX,0.,(Float_t) fNPixelsX);
602 fHis->AddAt(new TH1F(spdName.Data(),"SPD maps",
603 fNPixelsX,0.,(Float_t) fNPixelsX), i);
607 //____________________________________________
609 void AliITSsimulationSPDdubna::ResetHistograms()
612 // Reset histograms for this detector
615 for ( int i=0;i<fNPixelsZ;i++ ) {
616 //PH if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
617 if (fHis->At(i)) ((TH1F*)fHis->At(i))->Reset();