ITS/AliITSsimulationSPD.cxx

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