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

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


#include "AliRun.h"
#include "AliITS.h"
#include "AliITShit.h"
#include "AliITSdigit.h"
#include "AliITSmodule.h"
#include "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;
  fMapA2=0;
  fHis = 0;
  fNoise=0.;
  fBaseline=0.;
  fNPixelsZ=0;
  fNPixelsX=0;
}


//_____________________________________________________________________________

AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) {
  // standard constructor

      fHis = 0;
      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 yPrev = 1e+6;   

    Float_t zPitch = fSegmentation->Dpz(0);
    Float_t xPitch = fSegmentation->Dpx(0);
  
    TObjArray *fHits = mod->GetHits();
    Int_t nhits = fHits->GetEntriesFast();
    if (!nhits) return;

    //cout<<"len,wid,dy,nx,nz,pitchx,pitchz ="<<spdLength<<","<<spdWidth<<","<<fSegmentation->Dy()<<","<<fNPixelsX<<","<<fNPixelsZ<<","<<xPitch<<","<<zPitch<<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;
    //cout<<"SPD: module,nhits ="<<module<<","<<nhits<<endl;
    for (hit=0;hit<nhits;hit++) {
        AliITShit *iHit = (AliITShit*) fHits->At(hit);
	Int_t layer = iHit->GetLayer();
        Float_t yPix0 = -73; 
        if(layer == 1) yPix0 = -77; 

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

	//b.b.        if(iHit->StatusEntering()) idhit=hit;
        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();   // the momenta at the        
        Float_t py = iHit->GetPYL();   // each  GEANT step 
        Float_t pz = iHit->GetPZL();
        Float_t ptot = 1000*sqrt(px*px+py*py+pz*pz);
	*/

	Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
	                                           // vertex
        pmod *= 1000;


        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();      
	//cout<<"hit,status,y ="<<hit<<","<<status<<","<<yPix<<endl;      

	// Check boundaries
	if(zPix  > spdLength/2) {
	  //cout<<"!!!1 z outside ="<<zPix<<endl;
         zPix = spdLength/2 - 10;
	 //cout<<"!!!2 z outside ="<<zPix<<endl;
	}
	if(zPix  < 0 && zPix < -spdLength/2) {
	  //cout<<"!!!1 z outside ="<<zPix<<endl;
         zPix = -spdLength/2 + 10;
	 //cout<<"!!!2 z outside ="<<zPix<<endl;
	}
	if(xPix  > spdWidth/2) {
	  //cout<<"!!!1 x outside ="<<xPix<<endl;
         xPix = spdWidth/2 - 10;
	 //cout<<"!!!2 x outside ="<<xPix<<endl;
	}
	if(xPix  < 0 && xPix < -spdWidth/2) {
	  //cout<<"!!!1 x outside ="<<xPix<<endl;
         xPix = -spdWidth/2 + 10;
	 //cout<<"!!!2 x outside ="<<xPix<<endl;
	}
	Int_t trdown = 0;

	// enter Si or after event in Si
	if (status == 66 ) {  
           zPix0 = zPix;
           xPix0 = xPix;
           yPrev = yPix; 
	}   

	Float_t depEnergy = iHit->GetIonization();
	// skip if the input point to Si       

	if(depEnergy <= 0.) continue;        

	// if track returns to the opposite direction:
	if (yPix < 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 = TMath::Abs(yPix - yPrev);
	Float_t ydif0 = TMath::Abs(yPrev - yPix0);

	if(ydif < 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 (projDif < 5 ) {
	   drPath = (yPix-yPix0)*1.e-4;  
           drPath = TMath::Abs(drPath);        // drift path in cm
	   sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm        
	   sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
           nsteps = 1;
	}  

	if(projDif > 5) tang = ydif/projDif;
	Float_t dCharge = charge/nsteps;       // charge in e- for one step
	Float_t dZ = zdif/nsteps;
	Float_t dX = xdif/nsteps;

	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(projDif >= 5) {
	      Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
		drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm 
	      if(trdown == 0) {
		drPath = TMath::Abs(drPath) + ydif0*1.e-4;
	      }
	      if(trdown == 1) {
		drPath = ydif0*1.e-4 - TMath::Abs(drPath);
                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]) {
	     //b.b.	     tracks[j1]=-3;
	     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");

      fHis=new TObjArray(fNPixelsZ);
      TString spdName("spd_");
      for (Int_t i=0;i<fNPixelsZ;i++) {
	   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);
      }
}

//____________________________________________

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