ITS/AliITSsimulationSPD.cxx

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