ITS/AliITSsimulationSPD.cxx

   1 #include <iostream.h>
   2 #include <TRandom.h>
   3 #include <TH1.h>
   4 #include <TMath.h>
   5 #include <TString.h>
   6 #include <TParticle.h>
   7
   8
   9 #include "AliRun.h"
  10 #include "AliITS.h"
  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"
  18
  19
  20
  21
  22 ClassImp(AliITSsimulationSPD)
  23 ////////////////////////////////////////////////////////////////////////
  24 // Version: 0
  25 // Written by Boris Batyunya
  26 // December 20 1999
  27 //
  28 // AliITSsimulationSPD is the simulation of SPDs
  29 //________________________________________________________________________
  30
  31
  32 AliITSsimulationSPD::AliITSsimulationSPD()
  33 {
  34   // constructor
  35   fResponse = 0;
  36   fSegmentation = 0;
  37   fMapA2=0;
  38   fHis = 0;
  39   fNoise=0.;
  40   fBaseline=0.;
  41   fNPixelsZ=0;
  42   fNPixelsX=0;
  43 }
  44
  45
  46 //_____________________________________________________________________________
  47
  48 AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) {
  49   // standard constructor
  50
  51       fHis = 0;
  52       fResponse = resp;
  53       fSegmentation = seg;
  54
  55       fResponse->GetNoiseParam(fNoise,fBaseline);
  56
  57       fMapA2 = new AliITSMapA2(fSegmentation);
  58
  59       //
  60
  61       fNPixelsZ=fSegmentation->Npz();
  62       fNPixelsX=fSegmentation->Npx();
  63
  64 }
  65
  66 //_____________________________________________________________________________
  67
  68 AliITSsimulationSPD::~AliITSsimulationSPD() {
  69   // destructor
  70
  71   delete fMapA2;
  72
  73   if (fHis) {
  74      fHis->Delete();
  75      delete fHis;
  76   }
  77 }
  78
  79
  80 //__________________________________________________________________________
  81 AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
  82   //     Copy Constructor
  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;
  90   return;
  91 }
  92
  93 //_________________________________________________________________________
  94 AliITSsimulationSPD&
  95   AliITSsimulationSPD::operator=(const AliITSsimulationSPD &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;
 104   return *this;
 105   }
 106 //_____________________________________________________________________________
 107
 108 void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy)
 109 {
 110   // digitize module
 111
 112     const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
 113                                       // for 3.6 eV/pair
 114     const Float_t kconv = 10000.;     // cm -> microns
 115
 116     Float_t spdLength = fSegmentation->Dz();
 117     Float_t spdWidth = fSegmentation->Dx();
 118
 119     Float_t difCoef, dum;
 120     fResponse->DiffCoeff(difCoef,dum);
 121
 122     Float_t zPix0 = 1e+6;
 123     Float_t xPix0 = 1e+6;
 124     Float_t yPrev = 1e+6;
 125
 126     Float_t zPitch = fSegmentation->Dpz(0);
 127     Float_t xPitch = fSegmentation->Dpx(0);
 128
 129     TObjArray *fHits = mod->GetHits();
 130     Int_t nhits = fHits->GetEntriesFast();
 131     if (!nhits) return;
 132
 133     //cout<<"len,wid,dy,nx,nz,pitchx,pitchz ="<<spdLength<<","<<spdWidth<<","<<fSegmentation->Dy()<<","<<fNPixelsX<<","<<fNPixelsZ<<","<<xPitch<<","<<zPitch<<endl;
 134   //  Array of pointers to the label-signal list
 135
 136     Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;;
 137     Float_t  **pList = new Float_t* [maxNDigits];
 138     memset(pList,0,sizeof(Float_t*)*maxNDigits);
 139     Int_t indexRange[4] = {0,0,0,0};
 140
 141     // Fill detector maps with GEANT hits
 142     // loop over hits in the module
 143     static Bool_t first;
 144     Int_t lasttrack=-2;
 145     Int_t hit, iZi, jz, jx;
 146     //cout<<"SPD: module,nhits ="<<module<<","<<nhits<<endl;
 147     for (hit=0;hit<nhits;hit++) {
 148         AliITShit *iHit = (AliITShit*) fHits->At(hit);
 149         Int_t layer = iHit->GetLayer();
 150         Float_t yPix0 = -73;
 151         if(layer == 1) yPix0 = -77;
 152
 153         // work with the idtrack=entry number in the TreeH
 154         Int_t idhit,idtrack;
 155         mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);
 156         //Int_t idtrack=mod->GetHitTrackIndex(hit);
 157         // or store straight away the particle position in the array
 158         // of particles :
 159
 160         //b.b.        if(iHit->StatusEntering()) idhit=hit;
 161         Int_t itrack = iHit->GetTrack();
 162         Int_t dray = 0;
 163
 164         if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
 165
 166         //        Int_t parent = iHit->GetParticle()->GetFirstMother();
 167         Int_t partcode = iHit->GetParticle()->GetPdgCode();
 168
 169 //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
 170 //                      310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 171
 172         /*
 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);
 177         */
 178
 179         Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
 180                                                    // vertex
 181         pmod *= 1000;
 182
 183
 184         if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 185                                                  // at p < 6 MeV/c
 186
 187
 188         //  Get hit z and x(r*phi) cordinates for each module (detector)
 189         //  in local system.
 190
 191         Float_t zPix = kconv*iHit->GetZL();
 192         Float_t xPix = kconv*iHit->GetXL();
 193         Float_t yPix = kconv*iHit->GetYL();
 194
 195         // Get track status
 196         Int_t status = iHit->GetTrackStatus();
 197         //cout<<"hit,status,y ="<<hit<<","<<status<<","<<yPix<<endl;
 198
 199         // Check boundaries
 200         if(zPix  > spdLength/2) {
 201           //cout<<"!!!1 z outside ="<<zPix<<endl;
 202          zPix = spdLength/2 - 10;
 203          //cout<<"!!!2 z outside ="<<zPix<<endl;
 204         }
 205         if(zPix  < 0 && zPix < -spdLength/2) {
 206           //cout<<"!!!1 z outside ="<<zPix<<endl;
 207          zPix = -spdLength/2 + 10;
 208          //cout<<"!!!2 z outside ="<<zPix<<endl;
 209         }
 210         if(xPix  > spdWidth/2) {
 211           //cout<<"!!!1 x outside ="<<xPix<<endl;
 212          xPix = spdWidth/2 - 10;
 213          //cout<<"!!!2 x outside ="<<xPix<<endl;
 214         }
 215         if(xPix  < 0 && xPix < -spdWidth/2) {
 216           //cout<<"!!!1 x outside ="<<xPix<<endl;
 217          xPix = -spdWidth/2 + 10;
 218          //cout<<"!!!2 x outside ="<<xPix<<endl;
 219         }
 220         Int_t trdown = 0;
 221
 222         // enter Si or after event in Si
 223         if (status == 66 ) {
 224            zPix0 = zPix;
 225            xPix0 = xPix;
 226            yPrev = yPix;
 227         }
 228
 229         Float_t depEnergy = iHit->GetIonization();
 230         // skip if the input point to Si
 231
 232         if(depEnergy <= 0.) continue;
 233
 234         // if track returns to the opposite direction:
 235         if (yPix < yPrev) {
 236             trdown = 1;
 237         }
 238
 239
 240         // take into account the holes diffusion inside the Silicon
 241         // the straight line between the entrance and exit points in Si is
 242         // divided into the several steps; the diffusion is considered
 243         // for each end point of step and charge
 244         // is distributed between the pixels through the diffusion.
 245
 246
 247         //  ---------- the diffusion in Z (beam) direction -------
 248
 249         Float_t charge = depEnergy*kEnToEl;         // charge in e-
 250         Float_t drPath = 0.;
 251         Float_t tang = 0.;
 252         Float_t sigmaDif = 0.;
 253         Float_t zdif = zPix - zPix0;
 254         Float_t xdif = xPix - xPix0;
 255         Float_t ydif = TMath::Abs(yPix - yPrev);
 256         Float_t ydif0 = TMath::Abs(yPrev - yPix0);
 257
 258         if(ydif < 1) continue; // ydif is not zero
 259
 260         Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
 261
 262         Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
 263         Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
 264
 265         // number of the steps along the track:
 266         Int_t nsteps = ndZ;
 267         if(ndX > ndZ) nsteps = ndX;
 268         if(nsteps < 6) nsteps = 6;  // minimum number of the steps
 269
 270         if (projDif < 5 ) {
 271            drPath = (yPix-yPix0)*1.e-4;
 272            drPath = TMath::Abs(drPath);        // drift path in cm
 273            sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm
 274            sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
 275            nsteps = 1;
 276         }
 277
 278         if(projDif > 5) tang = ydif/projDif;
 279         Float_t dCharge = charge/nsteps;       // charge in e- for one step
 280         Float_t dZ = zdif/nsteps;
 281         Float_t dX = xdif/nsteps;
 282
 283         for (iZi = 1;iZi <= nsteps;iZi++) {
 284             Float_t dZn = iZi*dZ;
 285             Float_t dXn = iZi*dX;
 286             Float_t zPixn = zPix0 + dZn;
 287             Float_t xPixn = xPix0 + dXn;
 288
 289             if(projDif >= 5) {
 290               Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
 291                 drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
 292               if(trdown == 0) {
 293                 drPath = TMath::Abs(drPath) + ydif0*1.e-4;
 294               }
 295               if(trdown == 1) {
 296                 drPath = ydif0*1.e-4 - TMath::Abs(drPath);
 297                 drPath = TMath::Abs(drPath);
 298               }
 299               sigmaDif = difCoef*sqrt(drPath);
 300               sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
 301             }
 302
 303             zPixn = (zPixn + spdLength/2.);
 304             xPixn = (xPixn + spdWidth/2.);
 305             Int_t nZpix, nXpix;
 306             fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
 307             zPitch = fSegmentation->Dpz(nZpix);
 308             fSegmentation->GetPadTxz(xPixn,zPixn);
 309             // set the window for the integration
 310             Int_t jzmin = 1;
 311             Int_t jzmax = 3;
 312             if(nZpix == 1) jzmin =2;
 313             if(nZpix == fNPixelsZ) jzmax = 2;
 314
 315             Int_t jxmin = 1;
 316             Int_t jxmax = 3;
 317             if(nXpix == 1) jxmin =2;
 318             if(nXpix == fNPixelsX) jxmax = 2;
 319
 320             Float_t zpix = nZpix;
 321             Float_t dZright = zPitch*(zpix - zPixn);
 322             Float_t dZleft = zPitch - dZright;
 323
 324             Float_t xpix = nXpix;
 325             Float_t dXright = xPitch*(xpix - xPixn);
 326             Float_t dXleft = xPitch - dXright;
 327
 328             Float_t dZprev = 0.;
 329             Float_t dZnext = 0.;
 330             Float_t dXprev = 0.;
 331             Float_t dXnext = 0.;
 332
 333             for(jz=jzmin; jz <=jzmax; jz++) {
 334                 if(jz == 1) {
 335                   dZprev = -zPitch - dZleft;
 336                   dZnext = -dZleft;
 337                 }
 338                 if(jz == 2) {
 339                   dZprev = -dZleft;
 340                   dZnext = dZright;
 341                 }
 342                 if(jz == 3) {
 343                   dZprev = dZright;
 344                   dZnext = dZright + zPitch;
 345                 }
 346                 // kz changes from 1 to the fNofPixels(270)
 347                 Int_t kz = nZpix + jz -2;
 348
 349                 Float_t zArg1 = dZprev/sigmaDif;
 350                 Float_t zArg2 = dZnext/sigmaDif;
 351                 Float_t zProb1 = TMath::Erfc(zArg1);
 352                 Float_t zProb2 = TMath::Erfc(zArg2);
 353                 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
 354
 355
 356                 // ----------- holes diffusion in X(r*phi) direction  --------
 357
 358                 if(dZCharge > 1.) {
 359                   for(jx=jxmin; jx <=jxmax; jx++) {
 360                      if(jx == 1) {
 361                        dXprev = -xPitch - dXleft;
 362                        dXnext = -dXleft;
 363                      }
 364                      if(jx == 2) {
 365                        dXprev = -dXleft;
 366                        dXnext = dXright;
 367                      }
 368                      if(jx == 3) {
 369                        dXprev = dXright;
 370                        dXnext = dXright + xPitch;
 371                      }
 372                      Int_t kx = nXpix + jx -2;
 373
 374                      Float_t xArg1 = dXprev/sigmaDif;
 375                      Float_t xArg2 = dXnext/sigmaDif;
 376                      Float_t xProb1 = TMath::Erfc(xArg1);
 377                      Float_t xProb2 = TMath::Erfc(xArg2);
 378                      Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
 379
 380                      if(dXCharge > 1.) {
 381                        Int_t index = kz-1;
 382
 383                        if (first) {
 384                           indexRange[0]=indexRange[1]=index;
 385                           indexRange[2]=indexRange[3]=kx-1;
 386                           first=kFALSE;
 387                        }
 388
 389                        indexRange[0]=TMath::Min(indexRange[0],kz-1);
 390                        indexRange[1]=TMath::Max(indexRange[1],kz-1);
 391                        indexRange[2]=TMath::Min(indexRange[2],kx-1);
 392                        indexRange[3]=TMath::Max(indexRange[3],kx-1);
 393
 394                        // build the list of digits for this module
 395                        Double_t signal=fMapA2->GetSignal(index,kx-1);
 396                        signal+=dXCharge;
 397                        fMapA2->SetHit(index,kx-1,(double)signal);
 398                      }      // dXCharge > 1 e-
 399                   }       // jx loop
 400                 }       // dZCharge > 1 e-
 401             }        // jz loop
 402         }         // iZi loop
 403
 404         if (status == 65) {   // the step is inside of Si
 405            zPix0 = zPix;
 406            xPix0 = xPix;
 407         }
 408         yPrev = yPix;
 409
 410         if(dray == 0) {
 411             GetList(itrack,idhit,pList,indexRange);
 412         }
 413
 414         lasttrack=itrack;
 415     }   // hit loop inside the module
 416
 417
 418     // introduce the electronics effects and do zero-suppression
 419     ChargeToSignal(pList);
 420
 421     // clean memory
 422
 423     fMapA2->ClearMap();
 424
 425
 426 }
 427
 428 //---------------------------------------------
 429 void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
 430 {
 431   // lop over nonzero digits
 432
 433
 434   //set protection
 435   for(int k=0;k<4;k++) {
 436      if (indexRange[k] < 0) indexRange[k]=0;
 437   }
 438
 439   for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
 440     for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
 441
 442         Float_t signal=fMapA2->GetSignal(iz,ix);
 443
 444         if (!signal) continue;
 445
 446         Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
 447         if(!pList[globalIndex]){
 448
 449            //
 450            // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
 451            //
 452
 453            pList[globalIndex] = new Float_t [9];
 454
 455            // set list to -3
 456
 457            *pList[globalIndex] = -3.;
 458            *(pList[globalIndex]+1) = -3.;
 459            *(pList[globalIndex]+2) = -3.;
 460            *(pList[globalIndex]+3) =  0.;
 461            *(pList[globalIndex]+4) =  0.;
 462            *(pList[globalIndex]+5) =  0.;
 463            *(pList[globalIndex]+6) = -1.;
 464            *(pList[globalIndex]+7) = -1.;
 465            *(pList[globalIndex]+8) = -1.;
 466
 467
 468            *pList[globalIndex] = (float)label;
 469            *(pList[globalIndex]+3) = signal;
 470            *(pList[globalIndex]+6) = (float)idhit;
 471         }
 472         else{
 473
 474           // check the signal magnitude
 475
 476           Float_t highest = *(pList[globalIndex]+3);
 477           Float_t middle = *(pList[globalIndex]+4);
 478           Float_t lowest = *(pList[globalIndex]+5);
 479
 480           signal -= (highest+middle+lowest);
 481
 482           //
 483           //  compare the new signal with already existing list
 484           //
 485
 486           if(signal<lowest) continue; // neglect this track
 487
 488           if (signal>highest){
 489             *(pList[globalIndex]+5) = middle;
 490             *(pList[globalIndex]+4) = highest;
 491             *(pList[globalIndex]+3) = signal;
 492
 493             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
 494             *(pList[globalIndex]+1) = *pList[globalIndex];
 495             *pList[globalIndex] = label;
 496
 497             *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
 498             *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
 499             *(pList[globalIndex]+6) = idhit;
 500           }
 501           else if (signal>middle){
 502             *(pList[globalIndex]+5) = middle;
 503             *(pList[globalIndex]+4) = signal;
 504
 505             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
 506             *(pList[globalIndex]+1) = label;
 507
 508             *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
 509             *(pList[globalIndex]+7) = idhit;
 510           }
 511           else{
 512             *(pList[globalIndex]+5) = signal;
 513             *(pList[globalIndex]+2) = label;
 514             *(pList[globalIndex]+8) = idhit;
 515           }
 516         }
 517     } // end of loop pixels in x
 518   } // end of loop over pixels in z
 519
 520
 521 }
 522
 523
 524 //---------------------------------------------
 525 void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
 526 {
 527   // add noise and electronics, perform the zero suppression and add the
 528   // digit to the list
 529
 530   AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 531
 532
 533   TRandom *random = new TRandom();
 534   Float_t threshold = (float)fResponse->MinVal();
 535
 536   Int_t digits[3], tracks[3], hits[3],gi,j1;
 537   Float_t charges[3];
 538   Float_t electronics;
 539   Float_t signal,phys;
 540   for(Int_t iz=0;iz<fNPixelsZ;iz++){
 541     for(Int_t ix=0;ix<fNPixelsX;ix++){
 542       electronics = fBaseline + fNoise*random->Gaus();
 543       signal = (float)fMapA2->GetSignal(iz,ix);
 544       signal += electronics;
 545       gi =iz*fNPixelsX+ix; // global index
 546       if (signal > threshold) {
 547          digits[0]=iz;
 548          digits[1]=ix;
 549          digits[2]=1;
 550          for(j1=0;j1<3;j1++){
 551            if (pList[gi]) {
 552              //b.b.          tracks[j1]=-3;
 553              tracks[j1] = (Int_t)(*(pList[gi]+j1));
 554              hits[j1] = (Int_t)(*(pList[gi]+j1+6));
 555            }else {
 556              tracks[j1]=-2; //noise
 557              hits[j1] = -1;
 558            }
 559            charges[j1] = 0;
 560          }
 561
 562          if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
 563            tracks[1] = -3;
 564            hits[1] = -1;
 565            tracks[2] = -3;
 566            hits[2] = -1;
 567          }
 568          if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
 569            tracks[1] = -3;
 570            hits[1] = -1;
 571          }
 572          if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
 573            tracks[2] = -3;
 574            hits[2] = -1;
 575          }
 576          if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
 577            tracks[2] = -3;
 578            hits[2] = -1;
 579          }
 580
 581          phys=0;
 582          aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
 583       }
 584       if(pList[gi]) delete [] pList[gi];
 585     }
 586   }
 587   delete [] pList;
 588
 589 }
 590
 591
 592 //____________________________________________
 593
 594 void AliITSsimulationSPD::CreateHistograms()
 595 {
 596   // create 1D histograms for tests
 597
 598       printf("SPD - create histograms\n");
 599
 600       fHis=new TObjArray(fNPixelsZ);
 601       TString spdName("spd_");
 602       for (Int_t i=0;i<fNPixelsZ;i++) {
 603            Char_t pixelz[4];
 604            sprintf(pixelz,"%d",i+1);
 605            spdName.Append(pixelz);
 606            (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
 607                               fNPixelsX,0.,(Float_t) fNPixelsX);
 608       }
 609 }
 610
 611 //____________________________________________
 612
 613 void AliITSsimulationSPD::ResetHistograms()
 614 {
 615     //
 616     // Reset histograms for this detector
 617     //
 618
 619     for ( int i=0;i<fNPixelsZ;i++ ) {
 620         if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
 621     }
 622
 623 }
 624
 625
 626
 627
 628
 629
 630
 631
 632