ITS/AliITSsimulationSSD.cxx

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