[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;
    //    cout << "module, nhits ="<<module<<","<<nhits<< endl;

  //  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;
    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();
        Int_t dray = 0;
   
	if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE; 

	//        Int_t parent = iHit->GetParticle()->GetFirstMother();
        Int_t partcode = iHit->GetParticle()->GetPdgCode();

//  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
//                      310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda

        Float_t px = iHit->GetPXL();
        Float_t py = iHit->GetPYL();
        Float_t pz = iHit->GetPZL();
        Float_t pmod = 1000*sqrt(px*px+py*py+pz*pz);


        if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
                                                 // at p < 6 MeV/c


	//  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;  

	if(dray == 0) {
            GetList(itrack,idhit,pList,indexRange);
	}

	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;
	 }

	 if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
	   tracks[1] = -3;
           hits[1] = -1;
	   tracks[2] = -3;
           hits[2] = -1;
         } 
	 if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
	   tracks[1] = -3;
           hits[1] = -1;   
         } 
	 if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
	   tracks[2] = -3;
           hits[2] = -1;   
         } 
	 if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
	   tracks[2] = -3;
           hits[2] = -1;   
         } 

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