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