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

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


#include "AliRun.h"
#include "AliITS.h"
#include "AliITSMapA2.h" 
#include "AliITSsimulationSPD.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) {
  // standard 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 spdLength = fSegmentation->Dz();
    Float_t spdWidth = fSegmentation->Dx();

    Float_t difCoef, dum; 
    fResponse->DiffCoeff(difCoef,dum); 

    Float_t zPix0 = 1e+6;
    Float_t xPix0 = 1e+6;
    Float_t yPix0 = 1e+6;
    Float_t yPrev = 1e+6;   
    Float_t zP0 = 100.;
    Float_t xP0 = 100.;

    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 + fNPixelsX ;; 
    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
    static Bool_t first=kTRUE;
    Int_t lasttrack=-2;
    Int_t hit, iZi, jz, jx;
    for (hit=0;hit<nhits;hit++) {
        AliITShit *iHit = (AliITShit*) fHits->At(hit);
	Int_t layer = iHit->GetLayer();


	// work with the idtrack=entry number in the TreeH
	Int_t idhit,idtrack;
	mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);    
	//Int_t idtrack=mod->GetHitTrackIndex(hit);  
        // or store straight away the particle position in the array
	// of particles : 
        Int_t itrack = iHit->GetTrack();
   
	//  Get hit z and x(r*phi) cordinates for each module (detector)
	//  in local system.

	Float_t zPix = kconv*iHit->GetZL();
	Float_t xPix = kconv*iHit->GetXL();
	Float_t yPix = kconv*iHit->GetYL();

	// Get track status
	Int_t status = iHit->GetTrackStatus();      
      
	// Check boundaries
	if(TMath::Abs(zPix) > spdLength/2.) {
	  printf("!! Zpix outside = %f\n",zPix);
	  if(status == 66) zP0=100;
	  continue;
	} 


	if (TMath::Abs(xPix) > spdWidth/2.) {
	  printf("!! Xpix outside = %f\n",xPix);
	  if (status == 66) xP0=100;
	  continue;
	}  

	Float_t zP = (zPix + spdLength/2.)/1000.;  
	Float_t xP = (xPix + spdWidth/2.)/1000.;  

	Int_t 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;
             zP0 = zP;
             xP0 = xP;
	}
	// enter Si only
	if (layer == 2 && status == 66 && yPix < -71.) {    
             yPix0 = yPix;
             zP0 = zP;
             xP0 = xP;
	}       
	Float_t depEnergy = iHit->GetIonization();
	// skip if the input point to Si       
	if(depEnergy <= 0.) continue;        
	// skip if the input point is outside of Si, but the next
	// point is inside of Si
	if(zP0 > 90 || xP0 > 90) 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 -------

	Float_t charge = depEnergy*kEnToEl;         // charge in e-
	Float_t drPath = 0.;   
	Float_t tang = 0.;
	Float_t 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

	Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
	Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
	Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1; 

	// number of the steps along the track:
	Int_t 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;
	Float_t dCharge = charge/nsteps;       // charge in e- for one step
	Float_t dZ = zdif/nsteps;
	Float_t 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++) {
            Float_t dZn = iZi*dZ;
	    Float_t dXn = iZi*dX;
	    Float_t zPixn = zPix0 + dZn;
	    Float_t xPixn = xPix0 + dXn;

	    if(TMath::Abs(projDif) >= 5.) {
	      Float_t 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) {
		Float_t 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 + spdLength/2.);  
	    xPixn = (xPixn + spdWidth/2.);  
            Int_t nZpix, nXpix;
            fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
	    zPitch = fSegmentation->Dpz(nZpix);
            fSegmentation->GetPadTxz(xPixn,zPixn);
	    // set the window for the integration
	    Int_t jzmin = 1;  
	    Int_t jzmax = 3; 
	    if(nZpix == 1) jzmin =2;
	    if(nZpix == fNPixelsZ) jzmax = 2; 

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

	    Float_t zpix = nZpix; 
	    Float_t dZright = zPitch*(zpix - zPixn);
	    Float_t dZleft = zPitch - dZright;

	    Float_t xpix = nXpix; 
	    Float_t dXright = xPitch*(xpix - xPixn);
	    Float_t dXleft = xPitch - dXright;

	    Float_t dZprev = 0.;
	    Float_t dZnext = 0.;
	    Float_t dXprev = 0.;
	    Float_t 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;
		} 
		// kz changes from 1 to the fNofPixels(270)  
		Int_t kz = nZpix + jz -2; 

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


		// ----------- 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;
		     } 
		     Int_t kx = nXpix + jx -2;  

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

		     if(dXCharge > 1.) {
		       Int_t index = kz-1;

		       if (first) {
                          indexRange[0]=indexRange[1]=index;
                          indexRange[2]=indexRange[3]=kx-1;
                          first=kFALSE;
		       }

                       indexRange[0]=TMath::Min(indexRange[0],kz-1);
                       indexRange[1]=TMath::Max(indexRange[1],kz-1);
                       indexRange[2]=TMath::Min(indexRange[2],kx-1);
                       indexRange[3]=TMath::Max(indexRange[3],kx-1);

		       // build the list of digits for this module	
                       Double_t signal=fMapA2->GetSignal(index,kx-1);
                       signal+=dXCharge;
                       fMapA2->SetHit(index,kx-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 != itrack || hit==(nhits-1)) {
            GetList(itrack,idhit,pList,indexRange);
            first=kTRUE;
	}

	lasttrack=itrack;
    }   // 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,Int_t idhit,Float_t **pList,Int_t *indexRange)
{
  // lop over nonzero digits

   
  //set protection
  for(int k=0;k<4;k++) {
     if (indexRange[k] < 0) indexRange[k]=0;
  }

  for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
    for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){

        Float_t signal=fMapA2->GetSignal(iz,ix);
	if (!signal) continue;

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

           pList[globalIndex] = new Float_t [9];

	   // set list to -3 

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


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

	  // check the signal magnitude

          Float_t highest = *(pList[globalIndex]+3);
          Float_t middle = *(pList[globalIndex]+4);
          Float_t 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;

            *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
            *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
            *(pList[globalIndex]+6) = idhit;
	  }
          else if (signal>middle){
            *(pList[globalIndex]+5) = middle;
            *(pList[globalIndex]+4) = signal;

            *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
            *(pList[globalIndex]+1) = label;

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


}


//---------------------------------------------
void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
{
  // add noise and electronics, perform the zero suppression and add the
  // digit to the list

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

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

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

}


//____________________________________________

void AliITSsimulationSPD::CreateHistograms()
{
  // create 1D histograms for tests

      printf("SPD - create histograms\n");

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

}

//____________________________________________

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