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       fNPixelsZ=fSegmentation->Npz();
  61       fNPixelsX=fSegmentation->Npx();
  62
  63 }
  64
  65 //_____________________________________________________________________________
  66
  67 AliITSsimulationSPD::~AliITSsimulationSPD() {
  68   // destructor
  69
  70   delete fMapA2;
  71
  72   if (fHis) {
  73      fHis->Delete();
  74      delete fHis;
  75   }
  76 }
  77
  78
  79 //__________________________________________________________________________
  80 AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
  81   //     Copy Constructor
  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;
  89   return;
  90 }
  91
  92 //_________________________________________________________________________
  93 AliITSsimulationSPD&
  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;
 103   return *this;
 104   }
 105 //_____________________________________________________________________________
 106
 107 void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy)
 108 {
 109   // digitize module
 110
 111     const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
 112                                       // for 3.6 eV/pair
 113     const Float_t kconv = 10000.;     // cm -> microns
 114
 115     Float_t spdLength = fSegmentation->Dz();
 116     Float_t spdWidth = fSegmentation->Dx();
 117
 118     Float_t difCoef, dum;
 119     fResponse->DiffCoeff(difCoef,dum);
 120
 121     Float_t zPix0 = 1e+6;
 122     Float_t xPix0 = 1e+6;
 123     Float_t yPix0 = 1e+6;
 124     Float_t yPrev = 1e+6;
 125     Float_t zP0 = 100.;
 126     Float_t xP0 = 100.;
 127
 128     Float_t zPitch = fSegmentation->Dpz(0);
 129     Float_t xPitch = fSegmentation->Dpx(0);
 130
 131     //cout << "pitch per z: " << zPitch << endl;
 132     //cout << "pitch per r*phi: " << xPitch << endl;
 133
 134     TObjArray *fHits = mod->GetHits();
 135     Int_t nhits = fHits->GetEntriesFast();
 136     if (!nhits) return;
 137     //    cout << "module, nhits ="<<module<<","<<nhits<< endl;
 138
 139   //  Array of pointers to the label-signal list
 140
 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};
 145
 146     // Fill detector maps with GEANT hits
 147     // loop over hits in the module
 148     static Bool_t first;
 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();
 154
 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
 160         // of particles :
 161
 162         if(iHit->StatusEntering()) idhit=hit;
 163         Int_t itrack = iHit->GetTrack();
 164         Int_t dray = 0;
 165
 166         if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
 167
 168         //        Int_t parent = iHit->GetParticle()->GetFirstMother();
 169         Int_t partcode = iHit->GetParticle()->GetPdgCode();
 170
 171 //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
 172 //                      310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 173
 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);
 178
 179
 180         if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 181                                                  // at p < 6 MeV/c
 182
 183
 184         //  Get hit z and x(r*phi) cordinates for each module (detector)
 185         //  in local system.
 186
 187         Float_t zPix = kconv*iHit->GetZL();
 188         Float_t xPix = kconv*iHit->GetXL();
 189         Float_t yPix = kconv*iHit->GetYL();
 190
 191         // Get track status
 192         Int_t status = iHit->GetTrackStatus();
 193
 194         // Check boundaries
 195         if(TMath::Abs(zPix) > spdLength/2.) {
 196           printf("!! Zpix outside = %f\n",zPix);
 197           if(status == 66) zP0=100;
 198           continue;
 199         }
 200
 201
 202         if (TMath::Abs(xPix) > spdWidth/2.) {
 203           printf("!! Xpix outside = %f\n",xPix);
 204           if (status == 66) xP0=100;
 205           continue;
 206         }
 207
 208         Float_t zP = (zPix + spdLength/2.)/1000.;
 209         Float_t xP = (xPix + spdWidth/2.)/1000.;
 210
 211         Int_t trdown = 0;
 212
 213         // enter Si or after event in Si
 214         if (status == 66 ) {
 215            zPix0 = zPix;
 216            xPix0 = xPix;
 217            yPrev = yPix;
 218         }
 219         // enter Si only
 220         if (layer == 1 && status == 66 && yPix > 71.) {
 221              yPix0 = yPix;
 222              zP0 = zP;
 223              xP0 = xP;
 224         }
 225         // enter Si only
 226         if (layer == 2 && status == 66 && yPix < -71.) {
 227              yPix0 = yPix;
 228              zP0 = zP;
 229              xP0 = xP;
 230         }
 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) {
 239             yPix0 = yPrev;
 240             trdown = 1;
 241         }
 242         if (layer == 2 && yPix < yPrev) {
 243             yPix0 = yPrev;
 244             trdown = 1;
 245         }
 246
 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.
 252
 253
 254         //  ---------- the diffusion in Z (beam) direction -------
 255
 256         Float_t charge = depEnergy*kEnToEl;         // charge in e-
 257         Float_t drPath = 0.;
 258         Float_t tang = 0.;
 259         Float_t sigmaDif = 0.;
 260         Float_t zdif = zPix - zPix0;
 261         Float_t xdif = xPix - xPix0;
 262         Float_t ydif = yPix - yPix0;
 263
 264         if(TMath::Abs(ydif) < 0.1) continue; // Ydif is not zero
 265
 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;
 269
 270         // number of the steps along the track:
 271         Int_t nsteps = ndZ;
 272         if(ndX > ndZ) nsteps = ndX;
 273         if(nsteps < 6) nsteps = 6;  // minimum number of the steps
 274
 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;
 279
 280         if (TMath::Abs(projDif) < 5.0 ) {
 281            drPath = ydif*1.e-4;
 282            drPath = TMath::Abs(drPath);        // drift path in cm
 283            sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm
 284         }
 285
 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;
 291
 292             if(TMath::Abs(projDif) >= 5.) {
 293               Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
 294               if(trdown == 0) {
 295                 drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
 296                 drPath = TMath::Abs(drPath);
 297               }
 298               if(trdown == 1) {
 299                 Float_t dProjn = projDif/nsteps;
 300                 drPath = (projDif-(iZi-1)*dProjn)*tang*1.e-4;
 301                 drPath = TMath::Abs(drPath);
 302               }
 303               sigmaDif = difCoef*sqrt(drPath);
 304               sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
 305             }
 306             zPixn = (zPixn + spdLength/2.);
 307             xPixn = (xPixn + spdWidth/2.);
 308             Int_t nZpix, nXpix;
 309             fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
 310             zPitch = fSegmentation->Dpz(nZpix);
 311             fSegmentation->GetPadTxz(xPixn,zPixn);
 312             // set the window for the integration
 313             Int_t jzmin = 1;
 314             Int_t jzmax = 3;
 315             if(nZpix == 1) jzmin =2;
 316             if(nZpix == fNPixelsZ) jzmax = 2;
 317
 318             Int_t jxmin = 1;
 319             Int_t jxmax = 3;
 320             if(nXpix == 1) jxmin =2;
 321             if(nXpix == fNPixelsX) jxmax = 2;
 322
 323             Float_t zpix = nZpix;
 324             Float_t dZright = zPitch*(zpix - zPixn);
 325             Float_t dZleft = zPitch - dZright;
 326
 327             Float_t xpix = nXpix;
 328             Float_t dXright = xPitch*(xpix - xPixn);
 329             Float_t dXleft = xPitch - dXright;
 330
 331             Float_t dZprev = 0.;
 332             Float_t dZnext = 0.;
 333             Float_t dXprev = 0.;
 334             Float_t dXnext = 0.;
 335
 336             for(jz=jzmin; jz <=jzmax; jz++) {
 337                 if(jz == 1) {
 338                   dZprev = -zPitch - dZleft;
 339                   dZnext = -dZleft;
 340                 }
 341                 if(jz == 2) {
 342                   dZprev = -dZleft;
 343                   dZnext = dZright;
 344                 }
 345                 if(jz == 3) {
 346                   dZprev = dZright;
 347                   dZnext = dZright + zPitch;
 348                 }
 349                 // kz changes from 1 to the fNofPixels(270)
 350                 Int_t kz = nZpix + jz -2;
 351
 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;
 357
 358
 359                 // ----------- holes diffusion in X(r*phi) direction  --------
 360
 361                 if(dZCharge > 1.) {
 362                   for(jx=jxmin; jx <=jxmax; jx++) {
 363                      if(jx == 1) {
 364                        dXprev = -xPitch - dXleft;
 365                        dXnext = -dXleft;
 366                      }
 367                      if(jx == 2) {
 368                        dXprev = -dXleft;
 369                        dXnext = dXright;
 370                      }
 371                      if(jx == 3) {
 372                        dXprev = dXright;
 373                        dXnext = dXright + xPitch;
 374                      }
 375                      Int_t kx = nXpix + jx -2;
 376
 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;
 382
 383                      if(dXCharge > 1.) {
 384                        Int_t index = kz-1;
 385
 386                        if (first) {
 387                           indexRange[0]=indexRange[1]=index;
 388                           indexRange[2]=indexRange[3]=kx-1;
 389                           first=kFALSE;
 390                        }
 391
 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);
 396
 397                        // build the list of digits for this module
 398                        Double_t signal=fMapA2->GetSignal(index,kx-1);
 399                        signal+=dXCharge;
 400                        fMapA2->SetHit(index,kx-1,(double)signal);
 401                      }      // dXCharge > 1 e-
 402                   }       // jx loop
 403                 }       // dZCharge > 1 e-
 404             }        // jz loop
 405         }         // iZi loop
 406
 407         if (status == 65) {   // the step is inside of Si
 408            zPix0 = zPix;
 409            xPix0 = xPix;
 410         }
 411         yPrev = yPix;
 412
 413         if(dray == 0) {
 414             GetList(itrack,idhit,pList,indexRange);
 415         }
 416
 417         lasttrack=itrack;
 418     }   // hit loop inside the module
 419
 420
 421     // introduce the electronics effects and do zero-suppression
 422     ChargeToSignal(pList);
 423
 424     // clean memory
 425
 426     fMapA2->ClearMap();
 427
 428
 429 }
 430
 431 //---------------------------------------------
 432 void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
 433 {
 434   // lop over nonzero digits
 435
 436
 437   //set protection
 438   for(int k=0;k<4;k++) {
 439      if (indexRange[k] < 0) indexRange[k]=0;
 440   }
 441
 442   for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
 443     for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
 444
 445         Float_t signal=fMapA2->GetSignal(iz,ix);
 446
 447         if (!signal) continue;
 448
 449         Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
 450         if(!pList[globalIndex]){
 451
 452            //
 453            // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
 454            //
 455
 456            pList[globalIndex] = new Float_t [9];
 457
 458            // set list to -3
 459
 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.;
 469
 470
 471            *pList[globalIndex] = (float)label;
 472            *(pList[globalIndex]+3) = signal;
 473            *(pList[globalIndex]+6) = (float)idhit;
 474         }
 475         else{
 476
 477           // check the signal magnitude
 478
 479           Float_t highest = *(pList[globalIndex]+3);
 480           Float_t middle = *(pList[globalIndex]+4);
 481           Float_t lowest = *(pList[globalIndex]+5);
 482
 483           signal -= (highest+middle+lowest);
 484
 485           //
 486           //  compare the new signal with already existing list
 487           //
 488
 489           if(signal<lowest) continue; // neglect this track
 490
 491           if (signal>highest){
 492             *(pList[globalIndex]+5) = middle;
 493             *(pList[globalIndex]+4) = highest;
 494             *(pList[globalIndex]+3) = signal;
 495
 496             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
 497             *(pList[globalIndex]+1) = *pList[globalIndex];
 498             *pList[globalIndex] = label;
 499
 500             *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
 501             *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
 502             *(pList[globalIndex]+6) = idhit;
 503           }
 504           else if (signal>middle){
 505             *(pList[globalIndex]+5) = middle;
 506             *(pList[globalIndex]+4) = signal;
 507
 508             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
 509             *(pList[globalIndex]+1) = label;
 510
 511             *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
 512             *(pList[globalIndex]+7) = idhit;
 513           }
 514           else{
 515             *(pList[globalIndex]+5) = signal;
 516             *(pList[globalIndex]+2) = label;
 517             *(pList[globalIndex]+8) = idhit;
 518           }
 519         }
 520     } // end of loop pixels in x
 521   } // end of loop over pixels in z
 522
 523
 524 }
 525
 526
 527 //---------------------------------------------
 528 void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
 529 {
 530   // add noise and electronics, perform the zero suppression and add the
 531   // digit to the list
 532
 533   AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 534
 535
 536   TRandom *random = new TRandom();
 537   Float_t threshold = (float)fResponse->MinVal();
 538
 539   Int_t digits[3], tracks[3], hits[3],gi,j1;
 540   Float_t charges[3];
 541   Float_t electronics;
 542   Float_t signal,phys;
 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) {
 550          digits[0]=iz;
 551          digits[1]=ix;
 552          digits[2]=1;
 553          for(j1=0;j1<3;j1++){
 554            if (pList[gi]) {
 555              tracks[j1]=-3;
 556              tracks[j1] = (Int_t)(*(pList[gi]+j1));
 557              hits[j1] = (Int_t)(*(pList[gi]+j1+6));
 558            }else {
 559              tracks[j1]=-2; //noise
 560              hits[j1] = -1;
 561            }
 562            charges[j1] = 0;
 563          }
 564
 565          if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
 566            tracks[1] = -3;
 567            hits[1] = -1;
 568            tracks[2] = -3;
 569            hits[2] = -1;
 570          }
 571          if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
 572            tracks[1] = -3;
 573            hits[1] = -1;
 574          }
 575          if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
 576            tracks[2] = -3;
 577            hits[2] = -1;
 578          }
 579          if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
 580            tracks[2] = -3;
 581            hits[2] = -1;
 582          }
 583
 584          phys=0;
 585          aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
 586       }
 587       if(pList[gi]) delete [] pList[gi];
 588     }
 589   }
 590   delete [] pList;
 591
 592 }
 593
 594
 595 //____________________________________________
 596
 597 void AliITSsimulationSPD::CreateHistograms()
 598 {
 599   // create 1D histograms for tests
 600
 601       printf("SPD - create histograms\n");
 602
 603       fHis=new TObjArray(fNPixelsZ);
 604       TString spdName("spd_");
 605       for (Int_t i=0;i<fNPixelsZ;i++) {
 606            Char_t pixelz[4];
 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);
 611       }
 612 }
 613
 614 //____________________________________________
 615
 616 void AliITSsimulationSPD::ResetHistograms()
 617 {
 618     //
 619     // Reset histograms for this detector
 620     //
 621     for ( int i=0;i<fNPixelsZ;i++ ) {
 622         if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
 623     }
 624
 625 }
 626
 627
 628
 629
 630
 631
 632
 633
 634