[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPDbari.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 "AliITSsimulationSPDbari.h"
#include "AliITSsegmentation.h"
#include "AliITSresponse.h"


ClassImp(AliITSsimulationSPDbari)
////////////////////////////////////////////////////////////////////////
// Version: 0
// Written by Rocco Caliandro
// from a model developed with T. Virgili and R.A. Fini
// June 15 2000
//
// AliITSsimulationSPD is the simulation of SPDs
//
//________________________________________________________________________

AliITSsimulationSPDbari::AliITSsimulationSPDbari(){
  // constructor
  fResponse     = 0;
  fSegmentation = 0;
  fHis          = 0;
  fThresh       = 0.;
  fSigma        = 0.;
  fCouplCol     = 0.;
  fCouplRow     = 0.;
}
//_____________________________________________________________________________

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

      fResponse->Thresholds(fThresh,fSigma);
      fResponse->GetNoiseParam(fCouplCol,fCouplRow);
      
      fMapA2 = new AliITSMapA2(fSegmentation);
   
      // 
      fNPixelsZ=fSegmentation->Npz();
      fNPixelsX=fSegmentation->Npx();

}

//_____________________________________________________________________________

AliITSsimulationSPDbari::~AliITSsimulationSPDbari() { 
  // destructor

  delete fMapA2;

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

//__________________________________________________________________________
AliITSsimulationSPDbari::AliITSsimulationSPDbari(const AliITSsimulationSPDbari &source){
  //     Copy Constructor 
  if(&source == this) return;
  this->fMapA2 = source.fMapA2;
  this->fThresh = source.fThresh;
  this->fSigma = source.fSigma;
  this->fCouplCol = source.fCouplCol;
  this->fCouplRow = source.fCouplRow;
  this->fNPixelsX = source.fNPixelsX;
  this->fNPixelsZ = source.fNPixelsZ;
  this->fHis = source.fHis;
  return;
}

//_________________________________________________________________________
AliITSsimulationSPDbari& 
  AliITSsimulationSPDbari::operator=(const AliITSsimulationSPDbari &source) {
  //    Assignment operator
  if(&source == this) return *this;
  this->fMapA2 = source.fMapA2;
  this->fThresh = source.fThresh;
  this->fSigma = source.fSigma;
  this->fCouplCol = source.fCouplCol;
  this->fCouplRow = source.fCouplRow;
  this->fNPixelsX = source.fNPixelsX;
  this->fNPixelsZ = source.fNPixelsZ;
  this->fHis = source.fHis;
  return *this;
  }
//_____________________________________________________________________________

void AliITSsimulationSPDbari::DigitiseModule(AliITSmodule *mod, Int_t module,
                                             Int_t dummy) {
  // digitize module


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


  printf("Module %d (%d hits) \n",module+1,nhits);


  Int_t  number=10000;
  Int_t    *frowpixel = new Int_t[number];
  Int_t    *fcolpixel = new Int_t[number];
  Double_t *fenepixel = new Double_t[number];

  // Array of pointers to store the track index of the digits
  // leave +1, otherwise pList crashes when col=256, row=192 
    Int_t maxNdigits = fNPixelsX*fNPixelsZ+fNPixelsX+1;
  Float_t  **pList = new Float_t* [maxNdigits];
  memset(pList,0,sizeof(Float_t*)*maxNdigits);


  // noise setting
  SetFluctuations(pList);


  // loop over hits in the module
  Int_t hitpos;
      for (hitpos=0;hitpos<nhits;hitpos++) {  
 
     HitToDigit(mod,hitpos,module,frowpixel,fcolpixel,fenepixel,pList);


  }// end loop over digits

     CreateDigit(nhits,module,pList);

  // clean memory
  delete[] frowpixel;
  delete[] fcolpixel;
  delete[] fenepixel;
  fMapA2->ClearMap();
  delete [] pList;

}
//_____________________________________________________________________________

void AliITSsimulationSPDbari::UpdateMap( Int_t row, Int_t col, Double_t ene) {
//
// updates the Map of signal, adding the energy  (ene) released by the current track
//
      Double_t signal; 
      signal = fMapA2->GetSignal(row,col);
      signal += ene;
      fMapA2->SetHit(row,col,signal);
                                         
 }
//_____________________________________________________________________________  
void AliITSsimulationSPDbari::HitToDigit(AliITSmodule *mod, Int_t hitpos, Int_t module, 
                                        Int_t *frowpixel, Int_t *fcolpixel,
					Double_t *fenepixel, Float_t **pList) {
//
//  Steering function to determine the digits associated to a given hit (hitpos)
//  The digits are created by charge sharing (ChargeSharing) and by
//  capacitive coupling (SetCoupling). At all the created digits is associated
//  the track number of the hit (ntrack)
//


   static Float_t x1l,y1l,z1l;
   Float_t x2l,y2l,z2l,etot;
   Int_t layer,r1,r2,c1,c2,row,col,npixel = 0;
   Int_t ntrack;
   Double_t ene;
   const Float_t kconv = 10000.;     // cm -> microns

   TObjArray *fHits = mod->GetHits();
   AliITShit *hit = (AliITShit*) fHits->At(hitpos);
   layer = hit->GetLayer();
   etot=hit->GetIonization();
   ntrack=hit->GetTrack();


        if (hit->GetTrackStatus()==66) {
	      hit->GetPositionL(x1l,y1l,z1l);
          // positions shifted and converted in microns 
          x1l = x1l*kconv + fSegmentation->Dx()/2.;
          z1l = z1l*kconv + fSegmentation->Dz()/2.;
        }
        else {
	      hit->GetPositionL(x2l,y2l,z2l);	      
          // positions  shifted and converted in microns
          x2l = x2l*kconv + fSegmentation->Dx()/2.;
          z2l = z2l*kconv + fSegmentation->Dz()/2.;


          // to account for 83750 effective sensitive area
          // introduced in geometry (Dz=83600)
          if (z1l>fSegmentation->Dz()) return;
          if (z2l>fSegmentation->Dz()) return;
          // to preserve for hit having x>12800
          if (x1l>fSegmentation->Dx()) return;
          if (x2l>fSegmentation->Dx()) return;

          //Get the col and row number starting from 1
          // the x direction is not inverted for the second layer!!!
	      fSegmentation->GetPadIxz(x1l, z1l, c1, r1); 
	      fSegmentation->GetPadIxz(x2l, z2l, c2, r2);

          // to account for unexpected equal entrance and 
          // exit coordinates
          if (x1l==x2l) x2l=x2l+x2l*0.000001;
          if (z1l==z2l) z2l=z2l+z2l*0.000001;

          // to account for tracks at the edge of the sensitive area
          // of SPDs
          if (x1l<0) return;
          if (x2l<0) return;
          if (z1l<0) return;
          if (z2l<0) return;
	      

	      if ((r1==r2) && (c1==c2)) 
	      {
             // no charge sharing
	         npixel = 1;		 
		     frowpixel[npixel-1] = r1;
		     fcolpixel[npixel-1] = c1;
  		     fenepixel[npixel-1] = etot;
          }
	      else {
             // charge sharing
	         ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
		         	       npixel,frowpixel,fcolpixel,fenepixel);

          }
                  

      	  for (Int_t npix=0;npix<npixel;npix++)
	      {
		   row = frowpixel[npix];
		   col = fcolpixel[npix];
		   ene = fenepixel[npix];
		   UpdateMap(row,col,ene);                   
		   GetList(ntrack,pList,row,col); 
		   // Starting capacitive coupling effect
		   SetCoupling(row,col,ntrack,pList); 
	      }
	    x1l=x2l;
	    y1l=y2l;
	    z1l=z2l;	     	     
        }
}

//_________________________________________________________________________

void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
                    Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
				    Int_t r2,Float_t etot,
				    Int_t &npixel,Int_t *frowpixel,
				    Int_t *fcolpixel,Double_t *fenepixel){
  //
  //  Take into account the geometrical charge sharing when the track
  //  crosses more than one pixel.
  //
  //Begin_Html
  /*
  <img src="picts/ITS/barimodel_2.gif">
  </pre>
  <br clear=left>
  <font size=+2 color=red>
  <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
  </font>
  <pre>
  */
  //End_Html


   Float_t xa,za,xb,zb,dx,dz,dtot,dm,refr,refm,refc;
   Float_t refn=0.;
   Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
   Int_t   dirx,dirz,rb,cb;


   Int_t flag,flagrow,flagcol;
  
   Double_t epar;


   npixel = 0;
   xa = x1l;
   za = z1l;
   dx = TMath::Abs(x1l-x2l);
   dz = TMath::Abs(z1l-z2l);
   dtot = TMath::Sqrt((dx*dx)+(dz*dz));   
   dm = (x2l - x1l) / (z2l - z1l);

   dirx = (Int_t) ((x2l - x1l) / dx);
   dirz = (Int_t) ((z2l - z1l) / dz);
   
   
   // calculate the x coordinate of  the pixel in the next column    
   // and the z coordinate of  the pixel in the next row    

   Float_t xpos, zpos;

   fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos); 

   Float_t xsize = fSegmentation->Dpx(0);
   Float_t zsize = fSegmentation->Dpz(r1);

   if (dirx == 1) refr = xpos+xsize/2.;
             else refr = xpos-xsize/2.;

   if (dirz == 1) refn = zpos+zsize/2.;
             else refn = zpos-zsize/2.;

   
   flag = 0;
   flagrow = 0;
   flagcol = 0;
   do
   {
       
      // calculate the x coordinate of the intersection with the pixel
      // in the next cell in row  direction

      refm = (refn - z1l)*dm + x1l;
   
      // calculate the z coordinate of the intersection with the pixel
      // in the next cell in column direction 

      refc = (refr - x1l)/dm + z1l;
      
      
      arefm = refm * dirx;
      arefr = refr * dirx;
      arefn = refn * dirz;
      arefc = refc * dirz;
            

      if ((arefm < arefr) && (arefn < arefc)){
        	 
         // the track goes in the pixel in the next cell in row direction
	     xb = refm;
	     zb = refn;
	     cb = c1;
	     rb = r1 + dirz;
	     azb = zb * dirz;
         az2l = z2l * dirz;
	     if (rb == r2) flagrow=1;
	     if (azb > az2l) {
	        zb = z2l;
	        xb = x2l;
	     }     

         // shift to the pixel in the next cell in row direction
         Float_t zsizeNext = fSegmentation->Dpz(rb);
	     refn += zsizeNext*dirz;

      }
      else {
         
         // the track goes in the pixel in the next cell in column direction
	     xb = refr;
	     zb = refc;
	     cb = c1 + dirx;
	     rb = r1;
	     axb = xb * dirx;
         ax2l = x2l * dirx;
         if (cb == c2) flagcol=1;
	     if (axb > ax2l) {
	        zb = z2l;
	        xb = x2l;
	     }

         // shift to the pixel in the next cell in column direction
         Float_t xsizeNext = fSegmentation->Dpx(cb);
	     refr += xsizeNext*dirx;
        
      }
      
      //calculate the energy lost in the crossed pixel      
      epar = TMath::Sqrt((xb-xa)*(xb-xa)+(zb-za)*(zb-za)); 
      epar = etot*(epar/dtot);

      //store row, column and energy lost in the crossed pixel


      frowpixel[npixel] = r1;
      fcolpixel[npixel] = c1;
      fenepixel[npixel] = epar;
      npixel++;
 
      // the exit point of the track is reached
      if (epar == 0) flag = 1;
      if ((r1 == r2) && (c1 == c2)) flag = 1;
      if (flag!=1) {
        r1 = rb;
        c1 = cb;
        xa = xb;
        za = zb;
      }
   
   } while (flag==0);

}
//___________________________________________________________________________
void AliITSsimulationSPDbari::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
                                          Float_t **pList) {
   //
   //  Take into account the coupling between adiacent pixels.
   //  The parameters probcol and probrow are the fractions of the
   //  signal in one pixel shared in the two adjacent pixels along
   //  the column and row direction, respectively.
   //
   //Begin_Html
   /*
   <img src="picts/ITS/barimodel_3.gif">
   </pre>
   <br clear=left>
   <font size=+2 color=red>
   <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
   </font>
   <pre>
   */
   //End_Html


   Int_t j1,j2,flag=0;
   Double_t pulse1,pulse2;
                              

   j1 = row;
   j2 = col;
  
   pulse1 = fMapA2->GetSignal(row,col);
   pulse2 = pulse1;

   for (Int_t isign=-1;isign<=1;isign+=2)
   {

// loop in row direction
      
      do
      {
         j1 += isign;
         pulse1 *= fCouplRow;                  
      
         if ((j1 < 0) || (j1 > fNPixelsZ-1) || (pulse1 < fThresh))
         { 
	       pulse1 = fMapA2->GetSignal(row,col);
	       j1 = row;
	       flag = 1;
         }
          else{                
		   UpdateMap(j1,col,pulse1);                   
		   GetList(ntrack,pList,j1,col); 
           flag = 0;
	     }
	 
      } while(flag == 0);          
      
      
// loop in column direction
      
      do
      {
         j2 += isign;
         pulse2 *= fCouplCol;                  
      
         if ((j2 < 0) || (j2 > (fNPixelsX-1)) || (pulse2 < fThresh))
         {                
	       pulse2 = fMapA2->GetSignal(row,col);
	       j2 = col;
	       flag = 1;
         }
          else{                
		   UpdateMap(row,j2,pulse2);                   
		   GetList(ntrack,pList,row,j2); 
           flag = 0;
	     }
	 
      } while(flag == 0);          
   
   }

}
//___________________________________________________________________________
void AliITSsimulationSPDbari::CreateDigit(Int_t nhits, Int_t module, Float_t
**pList) {                                   
  //
  // The pixels are fired if the energy deposited inside them is above
  // the threshold parameter ethr. Fired pixed are interpreted as digits
  // and stored in the file digitfilename.
  //

   AliITS *aliITS  = (AliITS*)gAlice->GetModule("ITS");   
 
   
   Int_t digits[3];
   Int_t tracks[3];
   Int_t hits[3];
   Float_t charges[3]; 
   Int_t gi,j1;
   

   if (nhits > 0) {
    
     for (Int_t r=1;r<=fNPixelsZ;r++) {
        for (Int_t c=1;c<=fNPixelsX;c++) {
   
           // check if the deposited energy in a pixel is above the threshold 
           Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
           if ( signal > fThresh) {
	          digits[0] = r-1;  // digits starts from 0
              digits[1] = c-1;  // digits starts from 0
              digits[2] = 1;  
	          gi =r*fNPixelsX+c; // global index
	          for(j1=0;j1<3;j1++){
                 tracks[j1] = (Int_t)(*(pList[gi]+j1));
                 charges[j1] = 0;
              }
              Float_t phys = 0;        
	          aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
	          if(pList[gi]) delete [] pList[gi];
	       }//endif of threshold condition
        }
     }// enddo on pixels
    }
    
}
//_____________________________________________________________________________

void AliITSsimulationSPDbari::GetList(Int_t label,Float_t **pList,
                                      Int_t row, Int_t col) {
  // loop over nonzero digits

  Int_t ix = col;
  Int_t iz = row;
  Int_t globalIndex;
  Float_t signal;
  Float_t highest,middle,lowest;

          
  signal=fMapA2->GetSignal(iz,ix);

  globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 1


  if(!pList[globalIndex])
  {
     // 
     // Create new list (6 elements - 3 signals and 3 tracks + total sig)
     //

     pList[globalIndex] = new Float_t [6];


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

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


	  // check the signal magnitude
      highest = *(pList[globalIndex]+3);
      middle  = *(pList[globalIndex]+4);
      lowest  = *(pList[globalIndex]+5);


      signal -= (highest+middle+lowest);


	  //
	  //  compare the new signal with already existing list
	  //
      if(signal<lowest) return; // neglect this track

      if (signal>highest)
      {
         *(pList[globalIndex]+5) = middle;
         *(pList[globalIndex]+4) = highest;
         *(pList[globalIndex]+3) = signal;
         *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
         *(pList[globalIndex]+1) = *pList[globalIndex];
         *(pList[globalIndex]) = label;
	  }
        else if (signal>middle)
      {
         *(pList[globalIndex]+5) = middle;
         *(pList[globalIndex]+4) = signal;
         *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
         *(pList[globalIndex]+1) = label;
	  }
        else
      {
         *(pList[globalIndex]+5) = signal;
         *(pList[globalIndex]+2) = label;
	  }
  }    
}
//_________________________________________________________________________ 
void AliITSsimulationSPDbari::SetFluctuations(Float_t **pList) {
  //
  //  Set the electronic noise and threshold non-uniformities to all the
  //  pixels in a detector.
  //  The parameter fSigma is the squared sum of the sigma due to noise
  //  and the sigma of the threshold distribution among pixels.
  //
  //Begin_Html
  /*
  <img src="picts/ITS/barimodel_1.gif">
  </pre>
  <br clear=left>
  <font size=+2 color=red>
  <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
  </font>
  <pre>
  */
  //End_Html
  
  
  TRandom random; 
  Double_t signal;

  Int_t iz,ix;
  for(iz=1;iz<=fNPixelsZ;iz++){
    for(ix=1;ix<=fNPixelsX;ix++){
      signal = fSigma*random.Gaus(); 
      fMapA2->SetHit(iz,ix,signal);

      // insert in the label-signal list the pixels fired only by noise
      if ( signal > fThresh) {
        Int_t globalIndex = iz*fNPixelsX+ix; 
        pList[globalIndex] = new Float_t [6];
        *(pList[globalIndex]) = -2.;
        *(pList[globalIndex]+1) = -2.;
        *(pList[globalIndex]+2) = -2.;
        *(pList[globalIndex]+3) =  signal;
        *(pList[globalIndex]+4) =  0.;
        *(pList[globalIndex]+5) =  0.;
      }
    } // end of loop on pixels
  } // end of loop on pixels
  
 }
//____________________________________________

void AliITSsimulationSPDbari::CreateHistograms() {
  // CreateHistograms

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

}

//____________________________________________

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

}
Commit	Line	Data
0aedb77e	1	#include <iostream.h>
	2	#include <TRandom.h>
	3	#include <TH1.h>
	4	#include <TMath.h>
	5	#include <TString.h>
	6	#include <TParticle.h>
	7
	8
	9	#include "AliRun.h"
	10	#include "AliITS.h"
	11	#include "AliITShit.h"
	12	#include "AliITSdigit.h"
	13	#include "AliITSmodule.h"
	14	#include "AliITSMapA2.h"
	15	#include "AliITSsimulationSPDbari.h"
	16	#include "AliITSsegmentation.h"
	17	#include "AliITSresponse.h"
	18
	19
	20	ClassImp(AliITSsimulationSPDbari)
	21	////////////////////////////////////////////////////////////////////////
	22	// Version: 0
	23	// Written by Rocco Caliandro
	24	// from a model developed with T. Virgili and R.A. Fini
	25	// June 15 2000
	26	//
	27	// AliITSsimulationSPD is the simulation of SPDs
	28	//
	29	//________________________________________________________________________
	30
	31	AliITSsimulationSPDbari::AliITSsimulationSPDbari(){
	32	// constructor
	33	fResponse = 0;
	34	fSegmentation = 0;
	35	fHis = 0;
	36	fThresh = 0.;
	37	fSigma = 0.;
	38	fCouplCol = 0.;
	39	fCouplRow = 0.;
	40	}
	41	//_____________________________________________________________________________
	42
	43	AliITSsimulationSPDbari::AliITSsimulationSPDbari(AliITSsegmentation seg, AliITSresponse resp) {
	44	// constructor
	45	fResponse = resp;
	46	fSegmentation = seg;
	47
	48	fResponse->Thresholds(fThresh,fSigma);
	49	fResponse->GetNoiseParam(fCouplCol,fCouplRow);
	50
	51	fMapA2 = new AliITSMapA2(fSegmentation);
	52
	53	//
	54	fNPixelsZ=fSegmentation->Npz();
	55	fNPixelsX=fSegmentation->Npx();
	56
	57	}
	58
	59	//_____________________________________________________________________________
	60
	61	AliITSsimulationSPDbari::~AliITSsimulationSPDbari() {
	62	// destructor
	63
	64	delete fMapA2;
65
66	if (fHis) {
67	fHis->Delete();
68	delete fHis;
69	}
70	}
71
72	//__________________________________________________________________________
73	AliITSsimulationSPDbari::AliITSsimulationSPDbari(const AliITSsimulationSPDbari &source){
74	// Copy Constructor
75	if(&source == this) return;
76	this->fMapA2 = source.fMapA2;
77	this->fThresh = source.fThresh;
78	this->fSigma = source.fSigma;
79	this->fCouplCol = source.fCouplCol;
80	this->fCouplRow = source.fCouplRow;
81	this->fNPixelsX = source.fNPixelsX;
82	this->fNPixelsZ = source.fNPixelsZ;
83	this->fHis = source.fHis;
84	return;
85	}
86
87	//_________________________________________________________________________
88	AliITSsimulationSPDbari&
89	AliITSsimulationSPDbari::operator=(const AliITSsimulationSPDbari &source) {
90	// Assignment operator
91	if(&source == this) return *this;
92	this->fMapA2 = source.fMapA2;
93	this->fThresh = source.fThresh;
94	this->fSigma = source.fSigma;
95	this->fCouplCol = source.fCouplCol;
96	this->fCouplRow = source.fCouplRow;
97	this->fNPixelsX = source.fNPixelsX;
98	this->fNPixelsZ = source.fNPixelsZ;
99	this->fHis = source.fHis;
100	return *this;
101	}
102	//_____________________________________________________________________________
103
104	void AliITSsimulationSPDbari::DigitiseModule(AliITSmodule *mod, Int_t module,
105	Int_t dummy) {
106	// digitize module
107
108
109	TObjArray *fHits = mod->GetHits();
110	Int_t nhits = fHits->GetEntriesFast();
111	if (!nhits) return;
112
113
114	printf("Module %d (%d hits) \n",module+1,nhits);
115
116
117	Int_t number=10000;
118	Int_t *frowpixel = new Int_t[number];
119	Int_t *fcolpixel = new Int_t[number];
120	Double_t *fenepixel = new Double_t[number];
121
122	// Array of pointers to store the track index of the digits
123	// leave +1, otherwise pList crashes when col=256, row=192
124	Int_t maxNdigits = fNPixelsX*fNPixelsZ+fNPixelsX+1;
125	Float_t *pList = new Float_t [maxNdigits];
126	memset(pList,0,sizeof(Float_t)maxNdigits);
127
128
129	// noise setting
130	SetFluctuations(pList);
131
132
133	// loop over hits in the module
134	Int_t hitpos;
135	for (hitpos=0;hitpos<nhits;hitpos++) {
136
137	HitToDigit(mod,hitpos,module,frowpixel,fcolpixel,fenepixel,pList);
138
139
140	}// end loop over digits
141
142	CreateDigit(nhits,module,pList);
143
144	// clean memory
145	delete[] frowpixel;
146	delete[] fcolpixel;
147	delete[] fenepixel;
148	fMapA2->ClearMap();
149	delete [] pList;
150
151	}
152	//_____________________________________________________________________________
153
154	void AliITSsimulationSPDbari::UpdateMap( Int_t row, Int_t col, Double_t ene) {
155	//
156	// updates the Map of signal, adding the energy (ene) released by the current track
157	//
158	Double_t signal;
159	signal = fMapA2->GetSignal(row,col);
160	signal += ene;
161	fMapA2->SetHit(row,col,signal);
162
163	}
164	//_____________________________________________________________________________
165	void AliITSsimulationSPDbari::HitToDigit(AliITSmodule *mod, Int_t hitpos, Int_t module,
166	Int_t frowpixel, Int_t fcolpixel,
167	Double_t fenepixel, Float_t *pList) {
168	//
169	// Steering function to determine the digits associated to a given hit (hitpos)
170	// The digits are created by charge sharing (ChargeSharing) and by
171	// capacitive coupling (SetCoupling). At all the created digits is associated
172	// the track number of the hit (ntrack)
173	//
174
175
176	static Float_t x1l,y1l,z1l;
177	Float_t x2l,y2l,z2l,etot;
178	Int_t layer,r1,r2,c1,c2,row,col,npixel = 0;
179	Int_t ntrack;
180	Double_t ene;
181	const Float_t kconv = 10000.; // cm -> microns
182
183	TObjArray *fHits = mod->GetHits();
184	AliITShit hit = (AliITShit) fHits->At(hitpos);
185	layer = hit->GetLayer();
186	etot=hit->GetIonization();
187	ntrack=hit->GetTrack();
188
189
190	if (hit->GetTrackStatus()==66) {
191	hit->GetPositionL(x1l,y1l,z1l);
192	// positions shifted and converted in microns
193	x1l = x1l*kconv + fSegmentation->Dx()/2.;
194	z1l = z1l*kconv + fSegmentation->Dz()/2.;
195	}
196	else {
197	hit->GetPositionL(x2l,y2l,z2l);
198	// positions shifted and converted in microns
199	x2l = x2l*kconv + fSegmentation->Dx()/2.;
200	z2l = z2l*kconv + fSegmentation->Dz()/2.;
201
202
203	// to account for 83750 effective sensitive area
204	// introduced in geometry (Dz=83600)
205	if (z1l>fSegmentation->Dz()) return;
206	if (z2l>fSegmentation->Dz()) return;
207	// to preserve for hit having x>12800
208	if (x1l>fSegmentation->Dx()) return;
209	if (x2l>fSegmentation->Dx()) return;
210
211	//Get the col and row number starting from 1
212	// the x direction is not inverted for the second layer!!!
213	fSegmentation->GetPadIxz(x1l, z1l, c1, r1);
214	fSegmentation->GetPadIxz(x2l, z2l, c2, r2);
215
a82b8484	216	// to account for unexpected equal entrance and
	217	// exit coordinates
	218	if (x1l==x2l) x2l=x2l+x2l*0.000001;
	219	if (z1l==z2l) z2l=z2l+z2l*0.000001;
	220
	221	// to account for tracks at the edge of the sensitive area
	222	// of SPDs
	223	if (x1l<0) return;
	224	if (x2l<0) return;
	225	if (z1l<0) return;
	226	if (z2l<0) return;
0aedb77e	227
a82b8484	228
0aedb77e	229	if ((r1==r2) && (c1==c2))
	230	{
	231	// no charge sharing
	232	npixel = 1;
	233	frowpixel[npixel-1] = r1;
	234	fcolpixel[npixel-1] = c1;
	235	fenepixel[npixel-1] = etot;
	236	}
	237	else {
	238	// charge sharing
	239	ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
	240	npixel,frowpixel,fcolpixel,fenepixel);
	241
	242	}
	243
	244
	245	for (Int_t npix=0;npix<npixel;npix++)
	246	{
	247	row = frowpixel[npix];
	248	col = fcolpixel[npix];
	249	ene = fenepixel[npix];
	250	UpdateMap(row,col,ene);
	251	GetList(ntrack,pList,row,col);
	252	// Starting capacitive coupling effect
	253	SetCoupling(row,col,ntrack,pList);
	254	}
	255	x1l=x2l;
	256	y1l=y2l;
	257	z1l=z2l;
	258	}
	259	}
	260
	261	//_________________________________________________________________________
	262
	263	void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
	264	Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
	265	Int_t r2,Float_t etot,
	266	Int_t &npixel,Int_t *frowpixel,
	267	Int_t fcolpixel,Double_t fenepixel){
	268	//
	269	// Take into account the geometrical charge sharing when the track
	270	// crosses more than one pixel.
	271	//
	272	//Begin_Html
	273	/*
4b365332	274	<img src="picts/ITS/barimodel_2.gif">
0aedb77e	275	</pre>
	276	<br clear=left>
	277	<font size=+2 color=red>
	278	<a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
	279	</font>
	280	<pre>
	281	*/
	282	//End_Html
	283
	284
	285	Float_t xa,za,xb,zb,dx,dz,dtot,dm,refr,refm,refc;
	286	Float_t refn=0.;
	287	Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
	288	Int_t dirx,dirz,rb,cb;
	289
	290
	291	Int_t flag,flagrow,flagcol;
	292
	293	Double_t epar;
	294
	295
	296	npixel = 0;
	297	xa = x1l;
	298	za = z1l;
	299	dx = TMath::Abs(x1l-x2l);
	300	dz = TMath::Abs(z1l-z2l);
	301	dtot = TMath::Sqrt((dxdx)+(dzdz));
	302	dm = (x2l - x1l) / (z2l - z1l);
	303
	304	dirx = (Int_t) ((x2l - x1l) / dx);
	305	dirz = (Int_t) ((z2l - z1l) / dz);
	306
	307
	308	// calculate the x coordinate of the pixel in the next column
	309	// and the z coordinate of the pixel in the next row
	310
	311	Float_t xpos, zpos;
	312
	313	fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos);
	314
	315	Float_t xsize = fSegmentation->Dpx(0);
	316	Float_t zsize = fSegmentation->Dpz(r1);
	317
	318	if (dirx == 1) refr = xpos+xsize/2.;
	319	else refr = xpos-xsize/2.;
	320
	321	if (dirz == 1) refn = zpos+zsize/2.;
	322	else refn = zpos-zsize/2.;
	323
	324
	325	flag = 0;
	326	flagrow = 0;
	327	flagcol = 0;
	328	do
	329	{
	330
	331	// calculate the x coordinate of the intersection with the pixel
	332	// in the next cell in row direction
	333
	334	refm = (refn - z1l)*dm + x1l;
	335
	336	// calculate the z coordinate of the intersection with the pixel
	337	// in the next cell in column direction
	338
339	refc = (refr - x1l)/dm + z1l;
340
341
342	arefm = refm * dirx;
343	arefr = refr * dirx;
344	arefn = refn * dirz;
345	arefc = refc * dirz;
346
347
348	if ((arefm < arefr) && (arefn < arefc)){
349
350	// the track goes in the pixel in the next cell in row direction
351	xb = refm;
352	zb = refn;
353	cb = c1;
354	rb = r1 + dirz;
355	azb = zb * dirz;
356	az2l = z2l * dirz;
357	if (rb == r2) flagrow=1;
358	if (azb > az2l) {
359	zb = z2l;
360	xb = x2l;
361	}
362
363	// shift to the pixel in the next cell in row direction
364	Float_t zsizeNext = fSegmentation->Dpz(rb);
365	refn += zsizeNext*dirz;
366
367	}
368	else {
369
370	// the track goes in the pixel in the next cell in column direction
371	xb = refr;
372	zb = refc;
373	cb = c1 + dirx;
374	rb = r1;
375	axb = xb * dirx;
376	ax2l = x2l * dirx;
377	if (cb == c2) flagcol=1;
378	if (axb > ax2l) {
379	zb = z2l;
380	xb = x2l;
381	}
382
383	// shift to the pixel in the next cell in column direction
384	Float_t xsizeNext = fSegmentation->Dpx(cb);
385	refr += xsizeNext*dirx;
386
387	}
388
389	//calculate the energy lost in the crossed pixel
390	epar = TMath::Sqrt((xb-xa)(xb-xa)+(zb-za)(zb-za));
391	epar = etot*(epar/dtot);
392
393	//store row, column and energy lost in the crossed pixel
a82b8484	394
a82b8484	395
0aedb77e	396	frowpixel[npixel] = r1;
	397	fcolpixel[npixel] = c1;
	398	fenepixel[npixel] = epar;
	399	npixel++;
	400
	401	// the exit point of the track is reached
	402	if (epar == 0) flag = 1;
	403	if ((r1 == r2) && (c1 == c2)) flag = 1;
	404	if (flag!=1) {
	405	r1 = rb;
	406	c1 = cb;
	407	xa = xb;
	408	za = zb;
	409	}
	410
	411	} while (flag==0);
	412
	413	}
	414	//___________________________________________________________________________
	415	void AliITSsimulationSPDbari::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
	416	Float_t **pList) {
	417	//
	418	// Take into account the coupling between adiacent pixels.
	419	// The parameters probcol and probrow are the fractions of the
	420	// signal in one pixel shared in the two adjacent pixels along
	421	// the column and row direction, respectively.
	422	//
	423	//Begin_Html
	424	/*
4b365332	425	<img src="picts/ITS/barimodel_3.gif">
0aedb77e	426	</pre>
	427	<br clear=left>
	428	<font size=+2 color=red>
	429	<a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
	430	</font>
	431	<pre>
	432	*/
	433	//End_Html
	434
	435
	436	Int_t j1,j2,flag=0;
	437	Double_t pulse1,pulse2;
	438
	439
	440	j1 = row;
	441	j2 = col;
	442
	443	pulse1 = fMapA2->GetSignal(row,col);
	444	pulse2 = pulse1;
	445
	446	for (Int_t isign=-1;isign<=1;isign+=2)
	447	{
	448
	449	// loop in row direction
	450
	451	do
	452	{
	453	j1 += isign;
	454	pulse1 *= fCouplRow;
	455
	456	if ((j1 < 0) \|\| (j1 > fNPixelsZ-1) \|\| (pulse1 < fThresh))
	457	{
	458	pulse1 = fMapA2->GetSignal(row,col);
	459	j1 = row;
	460	flag = 1;
	461	}
	462	else{
	463	UpdateMap(j1,col,pulse1);
	464	GetList(ntrack,pList,j1,col);
	465	flag = 0;
	466	}
	467
	468	} while(flag == 0);
	469
	470
	471	// loop in column direction
	472
	473	do
	474	{
	475	j2 += isign;
	476	pulse2 *= fCouplCol;
	477
	478	if ((j2 < 0) \|\| (j2 > (fNPixelsX-1)) \|\| (pulse2 < fThresh))
	479	{
	480	pulse2 = fMapA2->GetSignal(row,col);
	481	j2 = col;
	482	flag = 1;
	483	}
	484	else{
	485	UpdateMap(row,j2,pulse2);
	486	GetList(ntrack,pList,row,j2);
	487	flag = 0;
	488	}
	489
490	} while(flag == 0);
491
492	}
493
494	}
495	//___________________________________________________________________________
496	void AliITSsimulationSPDbari::CreateDigit(Int_t nhits, Int_t module, Float_t
497	**pList) {
498	//
499	// The pixels are fired if the energy deposited inside them is above
500	// the threshold parameter ethr. Fired pixed are interpreted as digits
501	// and stored in the file digitfilename.
502	//
503
504	AliITS aliITS = (AliITS)gAlice->GetModule("ITS");
505
506
507	Int_t digits[3];
508	Int_t tracks[3];
509	Int_t hits[3];
510	Float_t charges[3];
511	Int_t gi,j1;
512
513
514	if (nhits > 0) {
515
516	for (Int_t r=1;r<=fNPixelsZ;r++) {
517	for (Int_t c=1;c<=fNPixelsX;c++) {
518
519	// check if the deposited energy in a pixel is above the threshold
520	Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
521	if ( signal > fThresh) {
522	digits[0] = r-1; // digits starts from 0
523	digits[1] = c-1; // digits starts from 0
524	digits[2] = 1;
525	gi =r*fNPixelsX+c; // global index
526	for(j1=0;j1<3;j1++){
527	tracks[j1] = (Int_t)(*(pList[gi]+j1));
528	charges[j1] = 0;
529	}
530	Float_t phys = 0;
531	aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
532	if(pList[gi]) delete [] pList[gi];
533	}//endif of threshold condition
534	}
535	}// enddo on pixels
536	}
537
538	}
539	//_____________________________________________________________________________
540
541	void AliITSsimulationSPDbari::GetList(Int_t label,Float_t **pList,
542	Int_t row, Int_t col) {
543	// loop over nonzero digits
544
545	Int_t ix = col;
546	Int_t iz = row;
547	Int_t globalIndex;
548	Float_t signal;
549	Float_t highest,middle,lowest;
550
551
552	signal=fMapA2->GetSignal(iz,ix);
553
554	globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 1
555
556
557	if(!pList[globalIndex])
558	{
559	//
560	// Create new list (6 elements - 3 signals and 3 tracks + total sig)
561	//
562
563	pList[globalIndex] = new Float_t [6];
564
565
566	// set list to -2
567	*(pList[globalIndex]) = -2.;
568	*(pList[globalIndex]+1) = -2.;
569	*(pList[globalIndex]+2) = -2.;
570	*(pList[globalIndex]+3) = 0.;
571	*(pList[globalIndex]+4) = 0.;
572	*(pList[globalIndex]+5) = 0.;
573
574	*pList[globalIndex] = (float)label;
575	*(pList[globalIndex]+3) = signal;
576	}
577	else{
578
579
580	// check the signal magnitude
581	highest = *(pList[globalIndex]+3);
582	middle = *(pList[globalIndex]+4);
583	lowest = *(pList[globalIndex]+5);
584
585
586	signal -= (highest+middle+lowest);
587
588
589	//
590	// compare the new signal with already existing list
591	//
592	if(signal<lowest) return; // neglect this track
593
594	if (signal>highest)
595	{
596	*(pList[globalIndex]+5) = middle;
597	*(pList[globalIndex]+4) = highest;
598	*(pList[globalIndex]+3) = signal;
599	(pList[globalIndex]+2) = (pList[globalIndex]+1);
600	(pList[globalIndex]+1) = pList[globalIndex];
601	*(pList[globalIndex]) = label;
602	}
603	else if (signal>middle)
604	{
605	*(pList[globalIndex]+5) = middle;
606	*(pList[globalIndex]+4) = signal;
607	(pList[globalIndex]+2) = (pList[globalIndex]+1);
608	*(pList[globalIndex]+1) = label;
609	}
610	else
611	{
612	*(pList[globalIndex]+5) = signal;
613	*(pList[globalIndex]+2) = label;
614	}
615	}
616	}
617	//_________________________________________________________________________
618	void AliITSsimulationSPDbari::SetFluctuations(Float_t **pList) {
619	//
620	// Set the electronic noise and threshold non-uniformities to all the
621	// pixels in a detector.
622	// The parameter fSigma is the squared sum of the sigma due to noise
623	// and the sigma of the threshold distribution among pixels.
624	//
625	//Begin_Html
626	/*
4b365332	627	<img src="picts/ITS/barimodel_1.gif">
0aedb77e	628	</pre>
	629	<br clear=left>
	630	<font size=+2 color=red>
	631	<a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
	632	</font>
	633	<pre>
	634	*/
	635	//End_Html
	636
	637
8813cd63	638	TRandom random;
0aedb77e	639	Double_t signal;
	640
	641	Int_t iz,ix;
	642	for(iz=1;iz<=fNPixelsZ;iz++){
	643	for(ix=1;ix<=fNPixelsX;ix++){
8813cd63	644	signal = fSigma*random.Gaus();
0aedb77e	645	fMapA2->SetHit(iz,ix,signal);
	646
	647	// insert in the label-signal list the pixels fired only by noise
	648	if ( signal > fThresh) {
	649	Int_t globalIndex = iz*fNPixelsX+ix;
	650	pList[globalIndex] = new Float_t [6];
	651	*(pList[globalIndex]) = -2.;
	652	*(pList[globalIndex]+1) = -2.;
	653	*(pList[globalIndex]+2) = -2.;
	654	*(pList[globalIndex]+3) = signal;
	655	*(pList[globalIndex]+4) = 0.;
	656	*(pList[globalIndex]+5) = 0.;
	657	}
	658	} // end of loop on pixels
	659	} // end of loop on pixels
	660
	661	}
	662	//____________________________________________
	663
	664	void AliITSsimulationSPDbari::CreateHistograms() {
	665	// CreateHistograms
	666
	667	Int_t i;
	668	for(i=0;i<fNPixelsZ;i++) {
8813cd63	669	TString spdname("spd_");
0aedb77e	670	Char_t candnum[4];
0aedb77e	671	sprintf(candnum,"%d",i+1);
8813cd63	672	spdname.Append(candnum);
8813cd63	673	(*fHis)[i] = new TH1F(spdname.Data(),"SPD maps",
0aedb77e	674	fNPixelsX,0.,(Float_t) fNPixelsX);
0aedb77e	675	}
	676
	677	}
	678
	679	//____________________________________________
	680
	681	void AliITSsimulationSPDbari::ResetHistograms() {
	682	//
	683	// Reset histograms for this detector
	684	//
	685	Int_t i;
	686	for(i=0;i<fNPixelsZ;i++ ) {
	687	if ((fHis)[i]) ((TH1F)(*fHis)[i])->Reset();
	688	}
	689
	690	}