[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx


#include <stdio.h>
#include <TObjArray.h>
#include <TRandom.h>

#include "AliITSsegmentationSSD.h"
#include "AliITSresponseSSD.h"
#include "AliITSsimulationSSD.h"
#include "AliITSdictSSD.h"
#include "AliITSdcsSSD.h"
#include "AliITS.h"
#include "AliITShit.h"
#include "AliITSdigit.h"
#include "AliITSmodule.h"
#include "AliRun.h"


ClassImp(AliITSsimulationSSD);
//------------------------------------------------------------
AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
                                         AliITSresponse *resp){
  // Constructor


    fSegmentation = seg;
    fResponse = resp;
    fDCS = new AliITSdcsSSD(seg,resp); 

    fNstrips = fSegmentation->Npx();
    fPitch = fSegmentation->Dpx(0);
    
    fP = new TArrayF(fNstrips+1); 
    fN = new TArrayF(fNstrips+1);
     
    fTracksP = new AliITSdictSSD[fNstrips+1];
    fTracksN = new AliITSdictSSD[fNstrips+1];

    
    fSteps  = 10;   // still hard-wired - set in SetDetParam and get it via  
                     // fDCS together with the others eventually    

}
//___________________________________________________________________________
AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD 
                                                                      &source){
// Operator =
    if(this==&source) return *this;

    this->fDCS = new AliITSdcsSSD(*(source.fDCS));
    this->fN   = new TArrayF(*(source.fN));
    this->fP   = new TArrayF(*(source.fP));
    this->fTracksP = new AliITSdictSSD(*(source.fTracksP));
    this->fTracksN = new AliITSdictSSD(*(source.fTracksN));
    this->fNstrips = source.fNstrips;
    this->fPitch   = source.fPitch;
    this->fSteps   = source.fSteps;
    return *this;
}
//_____________________________________________________________
AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
  // copy constructor
   *this = source;
}
//____________________________________________________________________________
AliITSsimulationSSD::~AliITSsimulationSSD() {
  // anihilator 


    if(fP) delete fP;
    if(fN) delete fN;

    if(fTracksP) delete [] fTracksP;
    if(fTracksN) delete [] fTracksN;

    delete fDCS;

} 
//_______________________________________________________________
//
// Hit to digit
//_______________________________________________________________
//
void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
                                         Int_t dummy) {
  // Digitizes one SSD module of hits.
    
    TObjArray *hits = mod->GetHits();
    Int_t nhits = hits->GetEntriesFast();
    if (!nhits) return;

    //printf("simSSD: module nhits %d %d\n",module,nhits);

    Int_t i;
    for(i=0; i<fNstrips; i++) {
       (*fP)[i] = 0;
       (*fN)[i] = 0;
       fTracksP[i].ZeroTracks();
       fTracksN[i].ZeroTracks();
    }
    
    for(i=0; i<nhits; i++) {
       Int_t idtrack=mod->GetHitTrackIndex(i);  
       HitToDigit(i,idtrack,nhits,hits);
    }
 
   
    ApplyNoise();
    ApplyCoupling();    
    ApplyThreshold();
    ApplyDAQ();
      

}

//---------------------------------------------------------------

void AliITSsimulationSSD::HitToDigit(Int_t & hitNo,Int_t idtrack,
                                     Int_t nhits,TObjArray *hits) {
  // Turns one or more hits in an SSD module into one or more digits.
     
    Int_t      stripP, stripN, i;
    Float_t    dsP, dsN;
    Float_t    sP, sN;
    Float_t    eP, eN;
    Float_t    arrayEP[10];         // hard-wired number of steps
    Float_t    arrayEN[10];
    Int_t      track = -3;
       
    Float_t    ionization = 0;
    Float_t    signal;
    
    AliITSdictSSD *dict;  
    
 
    // check if this is the right order !!!!!

    AliITShit *hitI = (AliITShit*)hits->At(hitNo++);
    AliITShit *hitE = (AliITShit*)hits->At(hitNo);


    while (!((hitE->StatusExiting()) || 
             (hitE->StatusDisappeared()) ||
             (hitE->StatusStop()))) {
       
        if (++hitNo<nhits) {
           ionization = hitE->GetIonization(); 
           hitE = (AliITShit*)hits->At(hitNo);
        }
    }   
    
        
    if (hitI->GetTrack() == hitE->GetTrack()) 
       //track = idtrack;
       track = hitI->GetTrack();
    else 
       printf("!!! Emergency !!!\n");
     
            
    ionization += hitE->GetIonization();
       
    const Float_t kconvm=10000.;  // cm -> microns

    Float_t xI, yI, zI;
    hitI->GetPositionL(xI, yI, zI);
    
    xI *= kconvm;
    yI *= kconvm;
    zI *= kconvm;
    
    Float_t xE, yE, zE;
    hitE->GetPositionL(xE, yE, zE);
    
    xE *= kconvm;
    yE *= kconvm;
    zE *= kconvm;

    Float_t dx = (xE - xI);
    Float_t dz = (zE - zI);
              
    
    // Debuging
    /*
    fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
   
       printf("%5d %8.3f %8.3f %8.3f %8.3f %d %d  %d\n", 
             hitNo, xI, zI, dx, dz, 
             stripP, stripN, track);
     printf("%10.5f %10d \n", ionization, hitI->fTrack); 
    */ 
    
    // end of debuging   
    
    
    eP=0;
    eN=0;
    
    for (i=0; i<fSteps; i++) {
        
      //        arrayEP[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
      //        arrayEN[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
        arrayEP[i] = ionization/fSteps;
        arrayEN[i] = ionization/fSteps;
       
        eP += arrayEP[i];
        eN += arrayEN[i];
    } 
       
    const Float_t kconv = 1.0e9 / 3.6;  // GeV -> e-hole pairs
       
    for(i=0; i<fSteps; i++) {
    
        arrayEP[i] = kconv * arrayEP[i] * (ionization / eP); 
        arrayEN[i] = kconv * arrayEN[i] * (ionization / eN);        
    }    
        
    dx /= fSteps;
    dz /= fSteps;  

    Float_t sigmaP, sigmaN; 
    fResponse->SigmaSpread(sigmaP,sigmaN); 

    //printf("SigmaP SigmaN %f %f\n",sigmaP, sigmaN);

    Float_t noiseP, noiseN;
    fResponse->GetNoiseParam(noiseP,noiseN);

    //printf("NoiseP NoiseN %f %f\n",noiseP, noiseN);

    for(i=0; i<fSteps; i++) {
        
        Int_t j;
       
        fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
        //printf("hitNo %d i xI zI stripP stripN %d %f %f %d %d\n",hitNo,i,xI, zI, stripP, stripN);
        dsP    = Get2Strip(1,stripP,xI, zI); // Between 0-1
        dsN    = Get2Strip(0,stripN,xI, zI); // Between 0-1

        //sP = sigmaP * sqrt(300. * i / (fSteps));
        //sN = sigmaN * sqrt(300. * i /(fSteps-i));

        sP = sigmaP * sqrt(300. * (i+1) / (fSteps));
        sN = sigmaN * sqrt(300. * (i+1) /(fSteps-i));


        sP = (i<2        && dsP>0.3 && dsP<0.7)? 20. : sP;  // square of (microns) 
        sN = (i>fSteps-2 && dsN>0.3 && dsN<0.7)? 20. : sN;  // square of (microns) 

        sP = (i==2 && dsP>0.4 && dsP<0.6)? 15. : sP;  // square of (microns) 
        sN = (i==8 && dsN>0.4 && dsN<0.6)? 15. : sN;  // square of (microns)        
        
        
	//printf("i=%d SigmaP SigmaN sP sN %f %f %e %e\n",i,sigmaP, sigmaN,sP,sN);
        
        for (j=-1; j<2; j++) {
            if (stripP+j<0 || stripP+j>fNstrips) continue;
            signal = arrayEP[i] * TMath::Abs( (F(j+0.5-dsP,sP)-F(j-0.5-dsP,sP)) );
            //printf("SimSSD::HitsToDigits:%d arrayEP[%d]=%e signal=%e\n",j,i,arrayEP[i],signal);
	    if (signal > noiseP/fSteps) {
               (*fP)[stripP+j] += signal;
               dict = (fTracksP+stripP+j);   
               (*dict).AddTrack(track);
            } 
	}  // end for j loop over neighboring strips

        for (j=-1; j<2; j++) {
            if (stripN+j<0 || stripN+j>fNstrips) continue;
            signal = arrayEN[i] * TMath::Abs( (F(j+0.5-dsN,sN)-F(j-0.5-dsN,sN)) );
            //printf("SimSSD::HitsToDigits:%d arrayEN[%d]=%e signal=%e\n",j,i,arrayEN[i],signal);
            if (signal > noiseN/fSteps) {
               (*fN)[stripN+j] += signal;
               dict = (fTracksN+stripN+j);    //co to jest
               (*dict).AddTrack(track);
            } 
        }  // end for j loop over neighboring strips
                
        xI += dx; 
        zI += dz; 
    }
    
    
}


//____________________________________________________________________
//
//  Private Methods for Simulation
//______________________________________________________________________
//

void AliITSsimulationSSD::ApplyNoise() {
  // Apply Noise.
   Float_t noiseP, noiseN;
   fResponse->GetNoiseParam(noiseP,noiseN);
       
    Int_t i;
    for(i = 0; i<fNstrips; i++) {
       (*fP)[i] += gRandom->Gaus(0,noiseP);
       (*fN)[i] += gRandom->Gaus(0,noiseN);
    }
}

//_________________________________________________________________________

void AliITSsimulationSSD::ApplyCoupling() {
  // Apply the effecto of electronic coupling between channels    
    Int_t i;
    for(i = 1; i<fNstrips-1; i++) {
      (*fP)[i] += (*fP)[i-1]*fDCS->GetCouplingPL() + (*fP)[i+1]*fDCS->GetCouplingPR();
      (*fN)[i] += (*fN)[i-1]*fDCS->GetCouplingNL() + (*fN)[i+1]*fDCS->GetCouplingNR();
    }
}

//__________________________________________________________________________

void AliITSsimulationSSD::ApplyThreshold() {
  // Applies the effect of a threshold on the signals for digitization.
   Float_t noiseP, noiseN;
   fResponse->GetNoiseParam(noiseP,noiseN);

   // or introduce the SetThresholds in fResponse  

    Int_t i;
    for(i=0; i<fNstrips; i++) {
       (*fP)[i] = ((*fP)[i] > noiseP*4) ? (*fP)[i] : 0;
       (*fN)[i] = ((*fN)[i] > noiseN*4) ? (*fN)[i] : 0; 
    }
        
}

//__________________________________________________________________________

void AliITSsimulationSSD::ApplyDAQ() {
  // Converts simulated signals to simulated ADC counts
    AliITS *its=(AliITS*)gAlice->GetModule("ITS");

    Float_t noiseP, noiseN;
    fResponse->GetNoiseParam(noiseP,noiseN);

    char opt[30],dummy[20];
    fResponse->ParamOptions(opt,dummy);

    Int_t i,j;
    if (strstr(opt,"SetInvalid")) {
      printf("invalid option %s\n",opt);
      // Set signal = 0 if invalid strip
      for(i=0; i<fNstrips; i++) {
         if (!(fDCS->IsValidP(i))) (*fP)[i] = 0;
	 if (!(fDCS->IsValidN(i))) (*fN)[i] = 0;
      }
    }
    
    Int_t digits[3], tracks[3], hits[3];
    Float_t charges[3];
    Float_t phys=0;
    for(i=0;i<3;i++) tracks[i]=-3;
    for(i=0; i<fNstrips; i++) { 
       if( (strstr(opt,"SetInvalid") && (*fP)[i] < noiseP*4) || !(*fP)[i]) continue;
          digits[0]=1;
          digits[1]=i;
          digits[2]=(int)(*fP)[i];
          for(j=0; j<(fTracksP+i)->GetNTracks(); j++) {
            if(j>2) continue;
	    if((fTracksP+i)->GetNTracks()) tracks[j]=(fTracksP+i)->GetTrack(j);
	    else tracks[j]=-2;
	    //printf("P side: i,j,tracks[j] %d %d %d\n",i,j,tracks[j]);
	    charges[j] = 0;
	    hits[j] = -1;
          }
          its->AddSimDigit(2,phys,digits,tracks,hits,charges);
          
          //cout << (fTracksP+i)->GetNTracks(); 
          //
	  //if ((fTracksP+i)->GetNTracks() == 0) {
          //   cout << d.fCoord2 << " " << d.fSignal << "\n"; 
          //}
    }
    
    
    for(i=0; i<fNstrips; i++) {
       if( (strstr(opt,"SetInvalid") && (*fN)[i] < noiseN*4)|| !(*fN)[i]) continue;
          digits[0]=0;
          digits[1]=i;
          digits[2]=(int)(*fN)[i];
          for( j=0; j<(fTracksN+i)->GetNTracks(); j++) {
            if(j>2) continue;
            if((fTracksN+i)->GetNTracks()) tracks[j]=(fTracksN+i)->GetTrack(j);
	    else tracks[j]=-2;
	    //printf("N side: i,j,tracks[j] %d %d %d\n",i,j,tracks[j]);
            charges[j] = 0;
	    hits[j] = -1;
          }
          its->AddSimDigit(2,phys,digits,tracks,hits,charges);
          
          //cout << (fTracksN+i)->GetNTracks();
          //if ((fTracksN+i)->GetNTracks() == 0) {
          //   cout << d.fCoord2 << " " << d.fSignal << "\n"; 
          //}
    }
    
}


//____________________________________________________________________________

Float_t AliITSsimulationSSD::F(Float_t x, Float_t s) {
  // Computes the integral of a gaussian at the mean valuse x with sigma s.
    //printf("SDD:F(%e,%e)\n",x,s);

    Float_t fval=0;
    if(s) fval=0.5*TMath::Erf(x * fPitch / s) ;
    else {
      Error("SSD simulation: F","sigma is zero!!!",s);
    }
    return fval;
} 

//______________________________________________________________________

Float_t AliITSsimulationSSD::Get2Strip(Int_t flag, Int_t iStrip, Float_t x, Float_t z){
  // Returns the relative space between two strips.

    // flag==1 for Pside, 0 for Nside

    Float_t stereoP, stereoN;
    fSegmentation->Angles(stereoP,stereoN);
    
    Float_t tanP=TMath::Tan(stereoP);
    Float_t tanN=TMath::Tan(stereoN);
 
    Float_t dx = fSegmentation->Dx();
    Float_t dz = fSegmentation->Dz();


     x += dx/2;
     z += dz/2; 
    
     if (flag) return (x - z*tanP) / fPitch - iStrip;       // from 0 to 1
     else  return (x - tanN*(dz - z)) / fPitch - iStrip;
}
//____________________________________________________________________________
Commit	Line	Data
b0f5e3fc	1
	2	#include <stdio.h>
	3	#include <TObjArray.h>
1ca7869b	4	#include <TRandom.h>
b0f5e3fc	5
9e8d6423	6	#include "AliITSsegmentationSSD.h"
9e8d6423	7	#include "AliITSresponseSSD.h"
b0f5e3fc	8	#include "AliITSsimulationSSD.h"
	9	#include "AliITSdictSSD.h"
	10	#include "AliITSdcsSSD.h"
	11	#include "AliITS.h"
9e8d6423	12	#include "AliITShit.h"
	13	#include "AliITSdigit.h"
	14	#include "AliITSmodule.h"
b0f5e3fc	15	#include "AliRun.h"
	16
	17
	18	ClassImp(AliITSsimulationSSD);
	19	//------------------------------------------------------------
	20	AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
	21	AliITSresponse *resp){
	22	// Constructor
	23
9e8d6423	24
b0f5e3fc	25	fSegmentation = seg;
	26	fResponse = resp;
	27	fDCS = new AliITSdcsSSD(seg,resp);
	28
	29	fNstrips = fSegmentation->Npx();
	30	fPitch = fSegmentation->Dpx(0);
	31
e8189707	32	fP = new TArrayF(fNstrips+1);
e8189707	33	fN = new TArrayF(fNstrips+1);
b0f5e3fc	34
e8189707	35	fTracksP = new AliITSdictSSD[fNstrips+1];
	36	fTracksN = new AliITSdictSSD[fNstrips+1];
	37
b0f5e3fc	38
	39	fSteps = 10; // still hard-wired - set in SetDetParam and get it via
	40	// fDCS together with the others eventually
	41
b0f5e3fc	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() {
9e8d6423	66	// anihilator
9e8d6423	67
e8189707	68
b0f5e3fc	69	if(fP) delete fP;
b0f5e3fc	70	if(fN) delete fN;
9e8d6423	71
	72	if(fTracksP) delete [] fTracksP;
	73	if(fTracksN) delete [] fTracksN;
b0f5e3fc	74
b0f5e3fc	75	delete fDCS;
9e8d6423	76
b0f5e3fc	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();
b0f5e3fc	89	if (!nhits) return;
b0f5e3fc	90
9e8d6423	91	//printf("simSSD: module nhits %d %d\n",module,nhits);
9e8d6423	92
b0f5e3fc	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
e8189707	107
b0f5e3fc	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;
e8189707	126	Float_t arrayEP[10]; // hard-wired number of steps
bba587ca	127	Float_t arrayEN[10];
9e8d6423	128	Int_t track = -3;
b0f5e3fc	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())
e8189707	154	//track = idtrack;
b0f5e3fc	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
b0f5e3fc	167	xI *= kconvm;
	168	yI *= kconvm;
	169	zI *= kconvm;
	170
	171	Float_t xE, yE, zE;
	172	hitE->GetPositionL(xE, yE, zE);
	173
b0f5e3fc	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	/*
e8189707	184	fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
b0f5e3fc	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;
b0f5e3fc	197
e8189707	198	for (i=0; i<fSteps; i++) {
b0f5e3fc	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
e8189707	234	fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
e8189707	235	//printf("hitNo %d i xI zI stripP stripN %d %f %f %d %d\n",hitNo,i,xI, zI, stripP, stripN);
b0f5e3fc	236	dsP = Get2Strip(1,stripP,xI, zI); // Between 0-1
	237	dsN = Get2Strip(0,stripN,xI, zI); // Between 0-1
	238
9e8d6423	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));
b0f5e3fc	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
9e8d6423	252
	253	//printf("i=%d SigmaP SigmaN sP sN %f %f %e %e\n",i,sigmaP, sigmaN,sP,sN);
	254
e8189707	255	for (j=-1; j<2; j++) {
b0f5e3fc	256	if (stripP+j<0 \|\| stripP+j>fNstrips) continue;
b0f5e3fc	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
e8189707	265
e8189707	266	for (j=-1; j<2; j++) {
b0f5e3fc	267	if (stripN+j<0 \|\| stripN+j>fNstrips) continue;
b0f5e3fc	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
b0f5e3fc	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] > noiseP4) ? (fP)[i] : 0;
	326	(fN)[i] = ((fN)[i] > noiseN4) ? (fN)[i] : 0;
b0f5e3fc	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
e8189707	340	char opt[30],dummy[20];
	341	fResponse->ParamOptions(opt,dummy);
	342
b0f5e3fc	343	Int_t i,j;
e8189707	344	if (strstr(opt,"SetInvalid")) {
9e8d6423	345	printf("invalid option %s\n",opt);
e8189707	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	}
b0f5e3fc	351	}
b0f5e3fc	352
e8189707	353	Int_t digits[3], tracks[3], hits[3];
b0f5e3fc	354	Float_t charges[3];
b0f5e3fc	355	Float_t phys=0;
9e8d6423	356	for(i=0;i<3;i++) tracks[i]=-3;
b0f5e3fc	357	for(i=0; i<fNstrips; i++) {
e8189707	358	if( (strstr(opt,"SetInvalid") && (fP)[i] < noiseP4) \|\| !(*fP)[i]) continue;
b0f5e3fc	359	digits[0]=1;
	360	digits[1]=i;
	361	digits[2]=(int)(*fP)[i];
	362	for(j=0; j<(fTracksP+i)->GetNTracks(); j++) {
e8189707	363	if(j>2) continue;
	364	if((fTracksP+i)->GetNTracks()) tracks[j]=(fTracksP+i)->GetTrack(j);
	365	else tracks[j]=-2;
9e8d6423	366	//printf("P side: i,j,tracks[j] %d %d %d\n",i,j,tracks[j]);
e8189707	367	charges[j] = 0;
e8189707	368	hits[j] = -1;
b0f5e3fc	369	}
e8189707	370	its->AddSimDigit(2,phys,digits,tracks,hits,charges);
b0f5e3fc	371
9e8d6423	372	//cout << (fTracksP+i)->GetNTracks();
	373	//
	374	//if ((fTracksP+i)->GetNTracks() == 0) {
	375	// cout << d.fCoord2 << " " << d.fSignal << "\n";
	376	//}
b0f5e3fc	377	}
	378
	379
	380	for(i=0; i<fNstrips; i++) {
e8189707	381	if( (strstr(opt,"SetInvalid") && (fN)[i] < noiseN4)\|\| !(*fN)[i]) continue;
b0f5e3fc	382	digits[0]=0;
	383	digits[1]=i;
	384	digits[2]=(int)(*fN)[i];
e8189707	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;
9e8d6423	389	//printf("N side: i,j,tracks[j] %d %d %d\n",i,j,tracks[j]);
e8189707	390	charges[j] = 0;
e8189707	391	hits[j] = -1;
b0f5e3fc	392	}
e8189707	393	its->AddSimDigit(2,phys,digits,tracks,hits,charges);
b0f5e3fc	394
9e8d6423	395	//cout << (fTracksN+i)->GetNTracks();
	396	//if ((fTracksN+i)->GetNTracks() == 0) {
	397	// cout << d.fCoord2 << " " << d.fSignal << "\n";
	398	//}
b0f5e3fc	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);
9e8d6423	409
	410	Float_t fval=0;
	411	if(s) fval=0.5TMath::Erf(x fPitch / s) ;
	412	else {
	413	Error("SSD simulation: F","sigma is zero!!!",s);
	414	}
	415	return fval;
b0f5e3fc	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	//____________________________________________________________________________