ITS/AliITSsimulationSSD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include <stdio.h>
  17 #include <TObjArray.h>
  18
  19 #include "AliITSsegmentationSSD.h"
  20 #include "AliITSresponseSSD.h"
  21 #include "AliITSsimulationSSD.h"
  22 #include "AliITSdictSSD.h"
  23 #include "AliITSdcsSSD.h"
  24 #include "AliITS.h"
  25 #include "AliRun.h"
  26
  27
  28 ClassImp(AliITSsimulationSSD);
  29 //------------------------------------------------------------
  30 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
  31                                          AliITSresponse *resp){
  32   // Constructor
  33
  34     fSegmentation = seg;
  35     fResponse = resp;
  36     fDCS = new AliITSdcsSSD(seg,resp);
  37
  38     fNstrips = fSegmentation->Npx();
  39     fPitch = fSegmentation->Dpx(0);
  40
  41     fP = new TArrayF(fNstrips);
  42     fN = new TArrayF(fNstrips);
  43
  44     fTracksP = new AliITSdictSSD[fNstrips];
  45     fTracksN = new AliITSdictSSD[fNstrips];
  46
  47     fSteps  = 10;   // still hard-wired - set in SetDetParam and get it via
  48                      // fDCS together with the others eventually
  49
  50
  51     //printf("SSD ctor: fNstrips fPitch %d %f\n",fNstrips, fPitch);
  52 }
  53 //___________________________________________________________________________
  54 AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD
  55                                                                       &source){
  56 // Operator =
  57     if(this==&source) return *this;
  58
  59     this->fDCS = new AliITSdcsSSD(*(source.fDCS));
  60     this->fN   = new TArrayF(*(source.fN));
  61     this->fP   = new TArrayF(*(source.fP));
  62     this->fTracksP = new AliITSdictSSD(*(source.fTracksP));
  63     this->fTracksN = new AliITSdictSSD(*(source.fTracksN));
  64     this->fNstrips = source.fNstrips;
  65     this->fPitch   = source.fPitch;
  66     this->fSteps   = source.fSteps;
  67     return *this;
  68 }
  69 //_____________________________________________________________
  70 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
  71   // copy constructor
  72    *this = source;
  73 }
  74 //____________________________________________________________________________
  75 AliITSsimulationSSD::~AliITSsimulationSSD() {
  76   // anihilator
  77
  78     if(fP) delete fP;
  79     if(fN) delete fN;
  80
  81     if(fTracksP) delete fTracksP;
  82     if(fTracksN) delete fTracksN;
  83
  84     delete fDCS;
  85 }
  86 //_______________________________________________________________
  87 //
  88 // Hit to digit
  89 //_______________________________________________________________
  90 //
  91 void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
  92                                          Int_t dummy) {
  93   // Digitizes one SSD module of hits.
  94
  95     TObjArray *hits = mod->GetHits();
  96     Int_t nhits = hits->GetEntriesFast();
  97     //printf("SSD: nhits %d\n",nhits);
  98     if (!nhits) return;
  99
 100     Int_t i;
 101     for(i=0; i<fNstrips; i++) {
 102        (*fP)[i] = 0;
 103        (*fN)[i] = 0;
 104        fTracksP[i].ZeroTracks();
 105        fTracksN[i].ZeroTracks();
 106     }
 107
 108     for(i=0; i<nhits; i++) {
 109        Int_t idtrack=mod->GetHitTrackIndex(i);
 110        HitToDigit(i,idtrack,nhits,hits);
 111     }
 112
 113
 114     ApplyNoise();
 115     ApplyCoupling();
 116     ApplyThreshold();
 117     ApplyDAQ();
 118
 119
 120 }
 121
 122 //---------------------------------------------------------------
 123
 124 void AliITSsimulationSSD::HitToDigit(Int_t & hitNo,Int_t idtrack,
 125                                      Int_t nhits,TObjArray *hits) {
 126   // Turns one or more hits in an SSD module into one or more digits.
 127
 128     Int_t      stripP, stripN, i;
 129     Float_t    dsP, dsN;
 130     Float_t    sP, sN;
 131     Float_t    eP, eN;
 132     Float_t    arrayEP[fSteps];
 133     Float_t    arrayEN[fSteps];
 134     Int_t      track = 0;
 135
 136     Float_t    ionization = 0;
 137     Float_t    signal;
 138
 139     AliITSdictSSD *dict;
 140
 141
 142     // check if this is the right order !!!!!
 143
 144     AliITShit *hitI = (AliITShit*)hits->At(hitNo++);
 145     AliITShit *hitE = (AliITShit*)hits->At(hitNo);
 146
 147
 148     while (!((hitE->StatusExiting()) ||
 149              (hitE->StatusDisappeared()) ||
 150              (hitE->StatusStop()))) {
 151
 152         if (++hitNo<nhits) {
 153            ionization = hitE->GetIonization();
 154            hitE = (AliITShit*)hits->At(hitNo);
 155         }
 156     }
 157
 158
 159     if (hitI->GetTrack() == hitE->GetTrack())
 160       //track = idtrack;
 161        track = hitI->GetTrack();
 162     else
 163        printf("!!! Emergency !!!\n");
 164
 165
 166     ionization += hitE->GetIonization();
 167
 168     const Float_t kconvm=10000.;  // cm -> microns
 169
 170     Float_t xI, yI, zI;
 171     hitI->GetPositionL(xI, yI, zI);
 172
 173     //Float_t zI =  hitI->GetZL();
 174
 175     xI *= kconvm;
 176     yI *= kconvm;
 177     zI *= kconvm;
 178
 179     Float_t xE, yE, zE;
 180     hitE->GetPositionL(xE, yE, zE);
 181
 182     //Float_t zE =  hitE->GetZL();
 183
 184     xE *= kconvm;
 185     yE *= kconvm;
 186     zE *= kconvm;
 187
 188     Float_t dx = (xE - xI);
 189     Float_t dz = (zE - zI);
 190
 191
 192     // Debuging
 193     /*
 194     fSegmentation->GetCellIxz(xI,zI,stripP,stripN);
 195
 196        printf("%5d %8.3f %8.3f %8.3f %8.3f %d %d  %d\n",
 197              hitNo, xI, zI, dx, dz,
 198              stripP, stripN, track);
 199      printf("%10.5f %10d \n", ionization, hitI->fTrack);
 200     */
 201
 202     // end of debuging
 203
 204
 205     eP=0;
 206     eN=0;
 207     //fNparticles++;
 208
 209     for(i=0; i<fSteps; i++) {
 210
 211       //        arrayEP[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
 212       //        arrayEN[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
 213         arrayEP[i] = ionization/fSteps;
 214         arrayEN[i] = ionization/fSteps;
 215
 216         eP += arrayEP[i];
 217         eN += arrayEN[i];
 218     }
 219
 220     const Float_t kconv = 1.0e9 / 3.6;  // GeV -> e-hole pairs
 221
 222     for(i=0; i<fSteps; i++) {
 223
 224         arrayEP[i] = kconv * arrayEP[i] * (ionization / eP);
 225         arrayEN[i] = kconv * arrayEN[i] * (ionization / eN);
 226     }
 227
 228     dx /= fSteps;
 229     dz /= fSteps;
 230
 231     Float_t sigmaP, sigmaN;
 232     fResponse->SigmaSpread(sigmaP,sigmaN);
 233
 234     //printf("SigmaP SigmaN %f %f\n",sigmaP, sigmaN);
 235
 236     Float_t noiseP, noiseN;
 237     fResponse->GetNoiseParam(noiseP,noiseN);
 238
 239     //printf("NoiseP NoiseN %f %f\n",noiseP, noiseN);
 240
 241     for(i=0; i<fSteps; i++) {
 242
 243         Int_t j;
 244
 245         fSegmentation->GetCellIxz(xI,zI,stripP,stripN);
 246         //printf("i xI zI stripP stripN %d %f %f %d %d\n",i,xI, zI, stripP, stripN);
 247         dsP    = Get2Strip(1,stripP,xI, zI); // Between 0-1
 248         dsN    = Get2Strip(0,stripN,xI, zI); // Between 0-1
 249
 250         sP = sigmaP * sqrt(300. * i / (fSteps));
 251         sN = sigmaN * sqrt(300. * i /(fSteps-i));
 252
 253
 254         sP = (i<2        && dsP>0.3 && dsP<0.7)? 20. : sP;  // square of (microns)
 255         sN = (i>fSteps-2 && dsN>0.3 && dsN<0.7)? 20. : sN;  // square of (microns)
 256
 257         sP = (i==2 && dsP>0.4 && dsP<0.6)? 15. : sP;  // square of (microns)
 258         sN = (i==8 && dsN>0.4 && dsN<0.6)? 15. : sN;  // square of (microns)
 259
 260
 261         //printf("i=%d SigmaP SigmaN sP sN %f %f %e %e\n",i,sigmaP, sigmaN,sP,sN);
 262
 263         for(j=-1; j<2; j++) {
 264
 265             if (stripP+j<0 || stripP+j>fNstrips) continue;
 266
 267             signal = arrayEP[i] * TMath::Abs( (F(j+0.5-dsP,sP)-F(j-0.5-dsP,sP)) );
 268             //printf("SimSSD::HitsToDigits:%d arrayEP[%d]=%e signal=%e\n",j,i,arrayEP[i],signal);
 269             if (signal > noiseP/fSteps) {
 270                (*fP)[stripP+j] += signal;
 271                dict = (fTracksP+stripP+j);
 272                (*dict).AddTrack(track);
 273             }
 274         }  // end for j loop over neighboring strips
 275         for(j=-1; j<2; j++) {
 276
 277             if (stripN+j<0 || stripN+j>fNstrips) continue;
 278
 279             signal = arrayEN[i] * TMath::Abs( (F(j+0.5-dsN,sN)-F(j-0.5-dsN,sN)) );
 280             //printf("SimSSD::HitsToDigits:%d arrayEN[%d]=%e signal=%e\n",j,i,arrayEN[i],signal);
 281             if (signal > noiseN/fSteps) {
 282                (*fN)[stripN+j] += signal;
 283                dict = (fTracksN+stripN+j);    //co to jest
 284                (*dict).AddTrack(track);
 285             }
 286         }  // end for j loop over neighboring strips
 287
 288         xI += dx;
 289         zI += dz;
 290     }
 291
 292
 293 }
 294
 295
 296 //____________________________________________________________________
 297 //
 298 //  Private Methods for Simulation
 299 //______________________________________________________________________
 300 //
 301
 302 void AliITSsimulationSSD::ApplyNoise() {
 303   // Apply Noise.
 304    Float_t noiseP, noiseN;
 305    fResponse->GetNoiseParam(noiseP,noiseN);
 306
 307     Int_t i;
 308     for(i = 0; i<fNstrips; i++) {
 309        (*fP)[i] += gRandom->Gaus(0,noiseP);
 310        (*fN)[i] += gRandom->Gaus(0,noiseN);
 311     }
 312 }
 313
 314 //_________________________________________________________________________
 315
 316 void AliITSsimulationSSD::ApplyCoupling() {
 317   // Apply the effecto of electronic coupling between channels
 318
 319     Int_t i;
 320     for(i = 1; i<fNstrips-1; i++) {
 321       (*fP)[i] += (*fP)[i-1]*fDCS->GetCouplingPL() + (*fP)[i+1]*fDCS->GetCouplingPR();
 322       (*fN)[i] += (*fN)[i-1]*fDCS->GetCouplingNL() + (*fN)[i+1]*fDCS->GetCouplingNR();
 323     }
 324 }
 325
 326 //__________________________________________________________________________
 327
 328 void AliITSsimulationSSD::ApplyThreshold() {
 329   // Applies the effect of a threshold on the signals for digitization.
 330    Float_t noiseP, noiseN;
 331    fResponse->GetNoiseParam(noiseP,noiseN);
 332
 333    // or introduce the SetThresholds in fResponse
 334
 335     Int_t i;
 336     for(i=0; i<fNstrips; i++) {
 337        (*fP)[i] = ((*fP)[i] > noiseP*4) ? (*fP)[i] : 0;
 338        (*fN)[i] = ((*fN)[i] > noiseN*4) ? (*fN)[i] : 0;
 339        //printf("SSD:(*fP)[i] (*fN)[i] %f %f \n",(*fP)[i], (*fN)[i]);
 340     }
 341
 342 }
 343
 344 //__________________________________________________________________________
 345
 346 void AliITSsimulationSSD::ApplyDAQ() {
 347   // Converts simulated signals to simulated ADC counts
 348     AliITS *its=(AliITS*)gAlice->GetModule("ITS");
 349
 350     Float_t noiseP, noiseN;
 351     fResponse->GetNoiseParam(noiseP,noiseN);
 352
 353     // Set signal = 0 if invalid strip
 354     Int_t i,j;
 355     for(i=0; i<fNstrips; i++) {
 356        if (!(fDCS->IsValidP(i))) (*fP)[i] = 0;
 357        if (!(fDCS->IsValidN(i))) (*fN)[i] = 0;
 358     }
 359
 360     Int_t digits[3], tracks[3];
 361     Float_t charges[3];
 362     Float_t phys=0;
 363     for(i=0; i<fNstrips; i++) {
 364        if ((*fP)[i] < noiseP*4) continue;
 365           digits[0]=1;
 366           digits[1]=i;
 367           digits[2]=(int)(*fP)[i];
 368           for(j=0; j<(fTracksP+i)->GetNTracks(); j++) {
 369              if(j>2) continue;
 370              tracks[j] = (fTracksP+i)->GetTrack(j);
 371              charges[j] = 0;
 372           }
 373           its->AddDigit(2,phys,digits,tracks,charges);
 374
 375           //cout << (fTracksP+i)->GetNTracks();
 376           //
 377           //if ((fTracksP+i)->GetNTracks() == 0) {
 378           //   cout << d.fCoord2 << " " << d.fSignal << "\n";
 379           //}
 380     }
 381
 382
 383     for(i=0; i<fNstrips; i++) {
 384        if ((*fN)[i] < noiseN*4) continue;
 385           digits[0]=0;
 386           digits[1]=i;
 387           digits[2]=(int)(*fN)[i];
 388           for(j=0; j<(fTracksN+i)->GetNTracks(); j++) {
 389              if(j>2) continue;
 390              tracks[j] = (fTracksN+i)->GetTrack(j);
 391              charges[j] = 0;
 392           }
 393           its->AddDigit(2,phys,digits,tracks,charges);
 394
 395           //cout << (fTracksN+i)->GetNTracks();
 396           //if ((fTracksN+i)->GetNTracks() == 0) {
 397           //   cout << d.fCoord2 << " " << d.fSignal << "\n";
 398           //}
 399     }
 400
 401 }
 402
 403
 404 //____________________________________________________________________________
 405
 406 Float_t AliITSsimulationSSD::F(Float_t x, Float_t s) {
 407   // Computes the integral of a gaussian at the mean valuse x with sigma s.
 408     //printf("SDD:F(%e,%e)\n",x,s);
 409     return 0.5*TMath::Erf(x * fPitch / s) ;
 410 }
 411
 412 //______________________________________________________________________
 413
 414 Float_t AliITSsimulationSSD::Get2Strip(Int_t flag, Int_t iStrip, Float_t x, Float_t z){
 415   // Returns the relative space between two strips.
 416
 417     // flag==1 for Pside, 0 for Nside
 418
 419     Float_t stereoP, stereoN;
 420     fSegmentation->Angles(stereoP,stereoN);
 421
 422     Float_t tanP=TMath::Tan(stereoP);
 423     Float_t tanN=TMath::Tan(stereoN);
 424
 425     Float_t dx = fSegmentation->Dx();
 426     Float_t dz = fSegmentation->Dz();
 427
 428
 429      x += dx/2;
 430      z += dz/2;
 431
 432      if (flag) return (x - z*tanP) / fPitch - iStrip;       // from 0 to 1
 433      else  return (x - tanN*(dz - z)) / fPitch - iStrip;
 434 }
 435 //____________________________________________________________________________