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

#include <iostream.h>
#include <TRandom.h>
#include <TH1.h>

#include "AliRun.h"
#include "AliITSMap.h"    // "AliITSMapA2.h"
#include "AliITSsimulationSPD.h"
#include "AliITSsegmentation.h"
#include "AliITSresponse.h"

ClassImp(AliITSsimulationSPD)
////////////////////////////////////////////////////////////////////////
// Version: 0
// Written by Boris Batyunya
// December 20 1999
//
// AliITSsimulationSPD is the simulation of SPDs
//
//________________________________________________________________________

AliITSsimulationSPD::AliITSsimulationSPD(){
  // constructor
  fResponse     = 0;
  fSegmentation = 0;
  fHis          = 0;
  fNoise        = 0.;
  fBaseline     = 0.;
}
//_____________________________________________________________________________

AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) {
  // constructor
      fResponse = resp;
      fSegmentation = seg;

      fResponse->GetNoiseParam(fNoise,fBaseline);

 
      fMapA2 = new AliITSMapA2(fSegmentation);
   
      // 
      fNPixelsZ=fSegmentation->Npz();
      fNPixelsX=fSegmentation->Npx();

}

//_____________________________________________________________________________

AliITSsimulationSPD::~AliITSsimulationSPD() { 
  // destructor

  delete fMapA2;

  if (fHis) {
     fHis->Delete(); 
     delete fHis;     
  }                
}

//__________________________________________________________________________
AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
  //     Copy Constructor 
  if(&source == this) return;
  this->fMapA2 = source.fMapA2;
  this->fNoise = source.fNoise;
  this->fBaseline = source.fBaseline;
  this->fNPixelsX = source.fNPixelsX;
  this->fNPixelsZ = source.fNPixelsZ;
  this->fHis = source.fHis;
  return;
}

//_________________________________________________________________________
AliITSsimulationSPD& 
  AliITSsimulationSPD::operator=(const AliITSsimulationSPD &source) {
  //    Assignment operator
  if(&source == this) return *this;
  this->fMapA2 = source.fMapA2;
  this->fNoise = source.fNoise;
  this->fBaseline = source.fBaseline;
  this->fNPixelsX = source.fNPixelsX;
  this->fNPixelsZ = source.fNPixelsZ;
  this->fHis = source.fHis;
  return *this;
  }
//_____________________________________________________________________________

void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy) {
  // digitize module
    const Float_t kEntoEl = 2.778e+8; // GeV->charge in electrons 
                                      // for 3.6 eV/pair 
    const Float_t kconv = 10000.;     // cm -> microns

    Float_t spdLenght = fSegmentation->Dz();
    Float_t spdWidth = fSegmentation->Dx();

    Float_t difCoef = fResponse->DiffCoeff();  

    Float_t zPix0 = 1e+6;
    Float_t xPix0 = 1e+6;
    Float_t yPix0 = 1e+6;
    Float_t yPrev = 1e+6;

    Float_t zPitch = fSegmentation->Dpz(0);
    Float_t xPitch = fSegmentation->Dpx(0);
  
    //cout << "pitch per z: " << zPitch << endl;
    //cout << "pitch per r*phi: " << xPitch << endl;

    TObjArray *fHits = mod->GetHits();
    Int_t nhits = fHits->GetEntriesFast();
    if (!nhits) return;


  //  Array of pointers to the label-signal list

    Int_t maxNdigits = fNPixelsX*fNPixelsZ; 
    Float_t  **pList = new Float_t* [maxNdigits]; 
    memset(pList,0,sizeof(Float_t*)*maxNdigits);
    Int_t indexRange[4] = {0,0,0,0};

    // Fill detector maps with GEANT hits
    // loop over hits in the module

    Int_t layer,idtrack,status,trdown,ndZ,ndX,nsteps,iZi,nZpix,nXpix;
    Int_t jzmin,jzmax,jxmin,jxmax,jx,jz,lz,lx,index;
    Float_t zArg1,zArg2,xArg1,xArg2;
    Float_t zProb1,zProb2,xProb1,xProb2;
    Float_t dZCharge,dXCharge;
    Double_t signal;
    Float_t projDif;  // RMS of xDif and zDif
    Float_t dProjn;  // RMS of dXn and dZn
    Float_t depEnergy; // The deposited energy from this hit
    Float_t zPix; // hit position in microns
    Float_t xPix; // hit position in microns
    Float_t yPix; // hit position in microns
    Float_t zPixn; // hit position in microns
    Float_t xPixn; // hit position in microns
    Float_t charge,dCharge;// charge in e-
    Float_t drPath;   
    Float_t tAng,dZ,dX,dXn,dZn;
    Float_t sigmaDif;
    Float_t dZright,dZleft;
    Float_t dXright,dXleft;
    Float_t dZprev,dZnext,dXprev,dXnext;

    static Bool_t first=kTRUE;
    Int_t lasttrack = -2,hit;
    for(hit=0;hit<nhits;hit++) {
        AliITShit *iHit = (AliITShit*) fHits->At(hit);
	layer = iHit->GetLayer();

	// work with the idtrack=entry number in the TreeH
	// Int_t idtrack=mod->GetHitTrackIndex(ii);  
        // or store straight away the particle position in the array
	// of particles : 
        idtrack = iHit->GetTrack();

	//if(module==11 || module==157) {
	  // Int_t track = iHit->fTrack;
	  // Int_t primary = gAlice->GetPrimary(track);
	  // Int_t parent = iHit->GetParticle()->GetFirstMother();
	  // printf("module,hit,track,primary,parent  %d %d %d %d %d \n",
	  //                                   md,hit,track,primary,parent);
	//}
   
	//  Get hit z and x(r*phi) cordinates for each module (detector)
	//  in local system.

	zPix = kconv*iHit->GetZL(); // Geant cm to microns
	xPix = kconv*iHit->GetXL(); // Geant cm to micron
        yPix = kconv*iHit->GetYL(); // Geant cm to micron

	// Get track status
	status = iHit->GetTrackStatus();      
	if(status != 66 && status != 68) {
	  printf("!!!!! no order status  %d\n",status);
	}

	trdown = 0;

	// enter Si or after event in Si
	if (status == 66 ) {  
           zPix0 = zPix;
           xPix0 = xPix;
           yPrev = yPix; 
	}   
	// enter Si only
	if (layer == 1 && status == 66 && yPix > 71.) {     
             yPix0 = yPix;
	}
	// enter Si only
	if (layer == 2 && status == 66 && yPix < -71.) {    
             yPix0 = yPix;
	}
	depEnergy = iHit->GetIonization();
	// skip if the input point to Si       
	if(depEnergy <= 0.) continue;  
	// if track returns to the opposite direction:
	if (layer == 1 && yPix > yPrev) {
	    yPix0 = yPrev;
            trdown = 1;
	}
	if (layer == 2 && yPix < yPrev) {
            yPix0 = yPrev;
            trdown = 1;
	}

	// take into account the holes diffusion inside the Silicon
	// the straight line between the entrance and exit points in Si is
	// divided into the several steps; the diffusion is considered 
	// for each end point of step and charge
	// is distributed between the pixels through the diffusion.
	

	//  ---------- the diffusion in Z (beam) direction -------

	charge = depEnergy*kEntoEl;         // charge in e-
	drPath = 0.;   
	tAng = 0.;
	sigmaDif = 0.; 
	Float_t zDif = zPix - zPix0;
	Float_t xDif = xPix - xPix0;
	Float_t yDif = yPix - yPix0;

	if(TMath::Abs(yDif) < 0.1) continue; // yDif is not zero


	projDif = sqrt(xDif*xDif + zDif*zDif);
	ndZ = (Int_t)TMath::Abs(zDif/zPitch) + 1;
	ndX = (Int_t)TMath::Abs(xDif/xPitch) + 1; 

	// number of the steps along the track:
	nsteps = ndZ;
	if(ndX > ndZ) nsteps = ndX;
	if(nsteps < 6) nsteps = 6;  // minimum number of the steps 

	if(TMath::Abs(projDif) > 5.0) tAng = yDif/projDif;
	dCharge = charge/nsteps;       // charge in e- for one step
	dZ = zDif/nsteps;
	dX = xDif/nsteps;

	if (TMath::Abs(projDif) < 5.0 ) {
	   drPath = yDif*1.e-4;  
           drPath = TMath::Abs(drPath);        // drift path in cm
	   sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm        
	}  

	for(iZi = 1;iZi <= nsteps;iZi++) {
            dZn = iZi*dZ;
	    dXn = iZi*dX;
	    zPixn = zPix0 + dZn;
	    xPixn = xPix0 + dXn;

	    if(TMath::Abs(projDif) >= 5.) {
	      dProjn = sqrt(dZn*dZn+dXn*dXn);
	      if(trdown == 0) {
		drPath = dProjn*tAng*1.e-4; // drift path for iZi step in cm 
		drPath = TMath::Abs(drPath);
	      }
	      if(trdown == 1) {
		dProjn = projDif/nsteps; 
		drPath = (projDif-(iZi-1)*dProjn)*tAng*1.e-4;
		drPath = TMath::Abs(drPath);
	      }
	      sigmaDif = difCoef*sqrt(drPath);    
	      sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
	    }
	    zPixn = (zPixn + spdLenght/2.);  
	    xPixn = (xPixn + spdWidth/2.);
            fSegmentation->GetCellIxz(xPixn,zPixn,nXpix,nZpix);
	    zPitch = fSegmentation->Dpz(nZpix);
	    // set the window for the integration
	    jzmin = 1;
	    jzmax = 3; 
	    if(nZpix == 1) jzmin =2;
	    if(nZpix == fNPixelsZ) jzmax = 2; 

	    jxmin = 1;  
	    jxmax = 3; 
	    if(nXpix == 1) jxmin =2;
	    if(nXpix == fNPixelsX) jxmax = 2; 

	    zPix    = nZpix;
	    dZright = zPitch*(zPix - zPixn);
	    dZleft  = zPitch - dZright;
	    xPix    = nXpix;
	    dXright = xPitch*(xPix - xPixn);
	    dXleft  = xPitch - dXright;
	    dZprev  = 0.;
	    dZnext  = 0.;
	    dXprev  = 0.;
	    dXnext  = 0.;

	    for(jz=jzmin; jz <=jzmax; jz++) {
	        if(jz == 1) {
		  dZprev = -zPitch - dZleft;
		  dZnext = -dZleft;
		} 
		if(jz == 2) {
		  dZprev = -dZleft;
		  dZnext = dZright;
		} 
		if(jz == 3) {
		  dZprev = dZright;
		  dZnext = dZright + zPitch;
		} 
		// lz changes from 1 to the fNofPixels(270)
		lz = nZpix + jz -2; 

		zArg1 = dZprev/sigmaDif;
		zArg2 = dZnext/sigmaDif;
		zProb1 = TMath::Erfc(zArg1);
		zProb2 = TMath::Erfc(zArg2);
		dZCharge =0.5*(zProb1-zProb2)*dCharge; 

		//printf("dZCharge %f \n",dZCharge);

		// ----------- holes diffusion in X(r*phi) direction  --------

		if(dZCharge > 1.) { 
		  for(jx=jxmin; jx <=jxmax; jx++) {
		     if(jx == 1) {
		       dXprev = -xPitch - dXleft;
		       dXnext = -dXleft;
		     } 
		     if(jx == 2) {
		       dXprev = -dXleft;
		       dXnext = dXright;
		     } 
		     if(jx == 3) {
		       dXprev = dXright;
		       dXnext = dXright + xPitch;
		     } 
		     lx = nXpix + jx -2;  

		     xArg1    = dXprev/sigmaDif;
		     xArg2    = dXnext/sigmaDif;
		     xProb1   = TMath::Erfc(xArg1);
		     xProb2   = TMath::Erfc(xArg2);
		     dXCharge =0.5*(xProb1-xProb2)*dZCharge; 

		     //printf("dXCharge %f \n",dXCharge);

		     if(dXCharge > 1.) {
		       index = lz-1;
		       if (first) {
                          indexRange[0]=indexRange[1]=index;
                          indexRange[2]=indexRange[3]=lx-1;
                          first=kFALSE;
		       }

                       indexRange[0]=TMath::Min(indexRange[0],lz-1);
                       indexRange[1]=TMath::Max(indexRange[1],lz-1);
                       indexRange[2]=TMath::Min(indexRange[2],lx-1);
                       indexRange[3]=TMath::Max(indexRange[3],lx-1);
		       // build the list of digits for this module	
                       signal=fMapA2->GetSignal(index,lx-1);
                       signal+=dXCharge;
                       fMapA2->SetHit(index,lx-1,(double)signal);
		     }      // dXCharge > 1 e-
		  }       // jx loop
		}       // dZCharge > 1 e-
	    }        // jz loop
	}         // iZi loop

        if (status == 65) {   // the step is inside of Si
	   zPix0 = zPix;
	   xPix0 = xPix;
        }
	yPrev = yPix;  //ch

	if (lasttrack != idtrack || hit==(nhits-1)) {
            GetList(idtrack,pList,indexRange);
            first=kTRUE;
	}
	lasttrack=idtrack;
    }   // hit loop inside the module

   
    // introduce the electronics effects and do zero-suppression
    ChargeToSignal(pList); 

    // clean memory

    fMapA2->ClearMap();


} 

//---------------------------------------------
void AliITSsimulationSPD::GetList(Int_t label,Float_t **pList,Int_t *indexRange) {
  // loop over nonzero digits

  Int_t ix,iz,globalIndex;
  Float_t signal;
  Float_t highest,middle,lowest;

  for(iz=indexRange[0];iz<indexRange[1]+1;iz++){
    for(ix=indexRange[2];ix<indexRange[3]+1;ix++){
      
      signal=fMapA2->GetSignal(iz,ix);

        globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 0!
        if(!pList[globalIndex]){
        
           // 
	   // Create new list (6 elements - 3 signals and 3 tracks + total sig)
	   //

           pList[globalIndex] = new Float_t [6];

	   // set list to -2 

	   *pList[globalIndex] = -2.;
	   *(pList[globalIndex]+1) = -2.;
	   *(pList[globalIndex]+2) = -2.;
	   *(pList[globalIndex]+3) =  0.;
	   *(pList[globalIndex]+4) =  0.;
	   *(pList[globalIndex]+5) =  0.;


	   *pList[globalIndex] = (float)label;
	   *(pList[globalIndex]+3) = signal;
        }
        else{

	  // check the signal magnitude

          highest = *(pList[globalIndex]+3);
          middle = *(pList[globalIndex]+4);
          lowest = *(pList[globalIndex]+5);

          signal -= (highest+middle+lowest);

	  //
	  //  compare the new signal with already existing list
	  //

          if(signal<lowest) continue; // neglect this track

          if (signal>highest){
            *(pList[globalIndex]+5) = middle;
            *(pList[globalIndex]+4) = highest;
            *(pList[globalIndex]+3) = signal;

            *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
            *(pList[globalIndex]+1) = *pList[globalIndex];
            *pList[globalIndex] = label;
	  }
          else if (signal>middle){
            *(pList[globalIndex]+5) = middle;
            *(pList[globalIndex]+4) = signal;

            *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
            *(pList[globalIndex]+1) = label;
	  }
          else{
            *(pList[globalIndex]+5) = signal;
            *(pList[globalIndex]+2) = label;
	  }
        }
    } // end of loop pixels in x
  } // end of loop over pixels in z


}


//---------------------------------------------
void AliITSsimulationSPD::ChargeToSignal(Float_t **pList) {
  // charge to signal  

  AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
  

  TRandom *random = new TRandom(); 
  Float_t threshold = (float)fResponse->MinVal();

  Int_t digits[3], tracks[3],gi,j1;
  Float_t charges[3];
  Float_t electronics;
  Float_t signal,phys;
  Int_t iz,ix;
  for(iz=0;iz<fNPixelsZ;iz++){
    for(ix=0;ix<fNPixelsX;ix++){
      electronics = fBaseline + fNoise*random->Gaus();
      signal = (float)fMapA2->GetSignal(iz,ix);
      signal += electronics;
      if (signal > threshold) {
	 digits[0]=iz;
	 digits[1]=ix;
	 digits[2]=1;
	 gi =iz*fNPixelsX+ix; // global index
	 for(j1=0;j1<3;j1++){
	   tracks[j1] = (Int_t)(*(pList[gi]+j1));
	   charges[j1] = 0;
	 }
         phys=0;
	 aliITS->AddDigit(0,phys,digits,tracks,charges);
         if(pList[gi]) delete [] pList[gi];
      }
    }
  }
  delete [] pList;


}


//____________________________________________

void AliITSsimulationSPD::CreateHistograms() {
  // CreateHistograms

      Int_t i;
      for(i=0;i<fNPixelsZ;i++) {
	   TString *spdname = new TString("spd_");
	   Char_t candnum[4];
	   sprintf(candnum,"%d",i+1);
	   spdname->Append(candnum);
	   (*fHis)[i] = new TH1F(spdname->Data(),"SPD maps",
                              fNPixelsX,0.,(Float_t) fNPixelsX);
	   delete spdname;
      }

}

//____________________________________________

void AliITSsimulationSPD::ResetHistograms() {
    //
    // Reset histograms for this detector
    //
    Int_t i;
    for(i=0;i<fNPixelsZ;i++ ) {
	if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
    }

}
Commit	Line	Data
b0f5e3fc	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(xDifxDif + zDifzDif);
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(dZndZn+dXndXn);
258	if(trdown == 0) {
259	drPath = dProjntAng1.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	}