Functions for bitio. Taken as is from The Data Compression Book

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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
Revision 1.13  2001/11/13 11:13:24  barbera
A protection against tracks with the same entrance and exit has been made more strict

Revision 1.12  2001/10/04 22:44:31  nilsen
Major changes in supppor of PreDigits (SDigits). Changes made with will make
it easier to suppor expected changes in AliITSHit class. Added use of new
class AliITSpList. Both SPD and SDD have added effects of Dead Channels. Both
of these will require addtional work as data bases of detectors and the like
are developed.

*/
#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 "AliITSpList.h"
#include "AliITSsimulationSPD.h"
#include "AliITSsegmentation.h"
#include "AliITSresponse.h"
#include "AliITSsegmentationSPD.h"
#include "AliITSresponseSPD.h"


ClassImp(AliITSsimulationSPD)
////////////////////////////////////////////////////////////////////////
// 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
//
//______________________________________________________________________
AliITSsimulationSPD::AliITSsimulationSPD(){
    // Default constructor

    fResponse     = 0;
    fSegmentation = 0;
    fHis          = 0;
    fMapA2        = 0;
/*
    fThresh       = 0.;
    fSigma        = 0.;
    fCouplCol     = 0.;
    fCouplRow     = 0.; */
}
//______________________________________________________________________
AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg,
					 AliITSresponse *resp) {
    // Standard constructor

    fResponse     = 0;
    fSegmentation = 0;
    fHis          = 0;
    fMapA2        = 0;
/*
    fThresh       = 0.;
    fSigma        = 0.;
    fCouplCol     = 0.;
    fCouplRow     = 0.*/
    Init((AliITSsegmentationSPD*)seg,(AliITSresponseSPD*)resp);
}
//______________________________________________________________________
void AliITSsimulationSPD::Init(AliITSsegmentationSPD *seg,
			       AliITSresponseSPD *resp) {
    // Initilizes the variables of AliITSsimulation SPD.

    fHis = 0;
    fResponse     = resp;
    fSegmentation = seg;
    fMapA2 = new AliITSMapA2(fSegmentation);
/*
    fResponse->Thresholds(fThresh,fSigma);
    fResponse->GetNoiseParam(fCouplCol,fCouplRow);
    fNPixelsZ = fSegmentation->Npz();
    fNPixelsX = fSegmentation->Npx();
*/
}
//______________________________________________________________________
AliITSsimulationSPD::~AliITSsimulationSPD() { 
    // destructor

    delete fMapA2;

    if (fHis) {
	fHis->Delete(); 
	delete fHis;     
    } // end if
}
//______________________________________________________________________
AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
    // Copy Constructor

    if(&source == this) return;

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

    if(&source == this) return *this;

    this->fMapA2    = source.fMapA2;
    this->fHis      = source.fHis;
/*
    this->fThresh   = source.fThresh;
    this->fSigma    = source.fSigma;
    this->fCouplCol = source.fCouplCol;
    this->fCouplRow = source.fCouplRow;
    this->fNPixelsX = source.fNPixelsX;
    this->fNPixelsZ = source.fNPixelsZ;
*/
    return *this;
}
//______________________________________________________________________
void AliITSsimulationSPD::SDigitiseModule(AliITSmodule *mod, Int_t dummy0,
                                             Int_t dummy1) {
    // Sum digitize module
    if (!(mod->GetNhits())) return; // if module has no hits then no Sdigits.
    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
    AliITSpList *pList = new AliITSpList(GetNPixelsZ()+1,GetNPixelsX()+1);

    HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,pList);

    WriteSDigits(pList);

    // clean memory
    delete[] frowpixel;
    delete[] fcolpixel;
    delete[] fenepixel;
    fMapA2->ClearMap();
    delete pList;
}
//______________________________________________________________________
void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t dummy0,
                                             Int_t dummy1) {
    // digitize module. Also need to digitize modules with only noise.

    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];
    Int_t    module     = mod->GetIndex();

    // Array of pointers to store the track index of the digits
    // leave +1, otherwise pList crashes when col=256, row=192 
    AliITSpList *pList = new AliITSpList(GetNPixelsZ()+1,GetNPixelsX()+1);

    // noise setting
    SetFluctuations(pList,module);

    HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,pList);

    // apply mask to SPD module
    SetMask();

    CreateDigit(module,pList);

    // clean memory
    delete[] frowpixel;
    delete[] fcolpixel;
    delete[] fenepixel;
    fMapA2->ClearMap();
    delete pList;
}
//______________________________________________________________________
void AliITSsimulationSPD::SDigitsToDigits(Int_t module,AliITSpList *pList) {
    // sum digits to Digits.

    FillMapFrompList(pList);

    // noise setting
    SetFluctuations(pList,module);

    // apply mask to SPD module
    SetMask();

    CreateDigit(module,pList);
}
//______________________________________________________________________
void AliITSsimulationSPD::UpdateMapSignal(Int_t row,Int_t col,Int_t trk,
					  Int_t hit,Int_t mod,Double_t ene,
					  AliITSpList *pList) {
    // updates the Map of signal, adding the energy  (ene) released by
    // the current track

    fMapA2->AddSignal(row,col,ene);
    pList->AddSignal(row,col,trk,hit,mod,ene);
}
//______________________________________________________________________
void AliITSsimulationSPD::UpdateMapNoise(Int_t row,Int_t col,Int_t mod,
					 Double_t ene,AliITSpList *pList) {
    // updates the Map of noise, adding the energy  (ene) give my noise

    fMapA2->AddSignal(row,col,ene);
    pList->AddNoise(row,col,mod,ene);
}
//______________________________________________________________________
void AliITSsimulationSPD::HitsToAnalogDigits(AliITSmodule *mod,
					     Int_t *frowpixel,Int_t *fcolpixel,
					     Double_t *fenepixel,
					     AliITSpList *pList) {
    // Loops over all hits to produce Analog/floting point digits. This
    // is also the first task in producing standard digits.
    
    // loop over hits in the module
    Int_t hitpos,nhits = mod->GetNhits();
    for (hitpos=0;hitpos<nhits;hitpos++) {
	HitToDigit(mod,hitpos,frowpixel,fcolpixel,fenepixel,pList);
    }// end loop over digits
}
//______________________________________________________________________
void AliITSsimulationSPD::HitToDigit(AliITSmodule *mod,Int_t hitpos,
				     Int_t *frowpixel,Int_t *fcolpixel,
				     Double_t *fenepixel,AliITSpList *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)
    Double_t x1l=0.0,y1l=0.0,z1l=0.0,x2l=0.0,y2l=0.0,z2l=0.0;
    Int_t r1,r2,c1,c2,row,col,npixel = 0;
    Int_t ntrack;
    Double_t ene=0.0,etot=0.0;
    const Float_t kconv = 10000.;     // cm -> microns
    const Float_t kconv1= 0.277e9;    // GeV -> electrons equivalent

    if(!(mod->LineSegmentL(hitpos,x1l,x2l,y1l,y2l,z1l,z2l,etot,ntrack)))return;

    x2l += x1l; y2l += y1l; z2l += z1l; // Convert to ending coordinate.
    // positions shifted and converted in microns
    x1l   = x1l*kconv + fSegmentation->Dx()/2.;
    z1l   = z1l*kconv + fSegmentation->Dz()/2.;
    // positions  shifted and converted in microns
    x2l   = x2l*kconv + fSegmentation->Dx()/2.;
    z2l   = z2l*kconv + fSegmentation->Dz()/2.;
    etot *= kconv1; // convert from GeV to electrons equivalent.
    Int_t module = mod->GetIndex();

    // to account for the effective sensitive area
    // introduced in geometry 
    if (z1l<0 || z1l>fSegmentation->Dz()) return;
    if (z2l<0 || z2l>fSegmentation->Dz()) return;
    if (x1l<0 || x1l>fSegmentation->Dx()) return;
    if (x2l<0 || 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.1;
    if (z1l==z2l) z2l=z2l+z2l*0.1;

    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);
    } // end if r1==r2 && c1==c2.

    for (Int_t npix=0;npix<npixel;npix++){
	row = frowpixel[npix];
	col = fcolpixel[npix];
	ene = fenepixel[npix];
	UpdateMapSignal(row,col,ntrack,hitpos,module,ene,pList); 
	// Starting capacitive coupling effect
	SetCoupling(row,col,ntrack,hitpos,module,pList); 
    } // end for npix
}
//______________________________________________________________________
void AliITSsimulationSPD::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);
    if (dx == 0.) dx = 0.01;
    dz     = TMath::Abs(z1l-z2l);
    if (dz == 0.) dz = 0.01;    
    dtot   = TMath::Sqrt((dx*dx)+(dz*dz));   
    dm     = (x2l - x1l) / (z2l - z1l);
    if (dm == 0.) dm = 0.01; 
    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-1);
    
    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;
	    } // end if
	    // shift to the pixel in the next cell in row direction
	    Float_t zsizeNext = fSegmentation->Dpz(rb-1);
	    //to account for cell at the borders of the detector
	    if(zsizeNext==0) zsizeNext = zsize;
	    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;
	    } // end ifaxb > ax2l

	    // shift to the pixel in the next cell in column direction
	    Float_t xsizeNext = fSegmentation->Dpx(cb-1);
	    //to account for cell at the borders of the detector
	    if(xsizeNext==0) xsizeNext = xsize;
	    refr += xsizeNext*dirx;
	} // end if (arefm < arefr) && (arefn < arefc)

	//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;
	} // end if flag!=1
    } while (flag==0);
}
//______________________________________________________________________
void AliITSsimulationSPD::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
				      Int_t idhit,Int_t module,
				      AliITSpList *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;
    Float_t couplR=0.0,couplC=0.0;

    GetCouplings(couplR,couplC);
    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 *= couplR;
	    if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
		pulse1 = fMapA2->GetSignal(row,col);
		j1 = row;
		flag = 1;
	    }else{
		UpdateMapSignal(j1,col,ntrack,idhit,module,pulse1,pList);
		flag = 0;
	    } // end if
	} while(flag == 0);
	// loop in column direction
      do{
	  j2 += isign;
	  pulse2 *= couplC;
	  if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
	      pulse2 = fMapA2->GetSignal(row,col);
	      j2 = col;
	      flag = 1;
	  }else{
	      UpdateMapSignal(row,j2,ntrack,idhit,module,pulse2,pList);
	      flag = 0;
	  } // end if
      } while(flag == 0);
    } // for isign
}
//______________________________________________________________________
void AliITSsimulationSPD::CreateDigit(Int_t module,AliITSpList *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. One also needs to write out
    // cases when there is only noise (nhits==0).

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

    Int_t digits[3];
    Int_t tracks[3];
    Int_t hits[3];
    Float_t charges[3]; 
    Int_t j1;

    for (Int_t r=1;r<=GetNPixelsZ();r++) {
	for (Int_t c=1;c<=GetNPixelsX();c++) {
	    // check if the deposited energy in a pixel is above the
	    // threshold 
	    Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
	    if ( signal > GetThreshold()) {
		digits[0] = r-1;  // digits starts from 0
		digits[1] = c-1;  // digits starts from 0
		//digits[2] = 1;  
		digits[2] =  (Int_t) signal;  // the signal is stored in
                                              //  electrons
		for(j1=0;j1<3;j1++){
		    tracks[j1] = pList->GetTrack(r,c,j1);
		    hits[j1]   = pList->GetHit(r,c,j1);
		    charges[j1] = 0;
		} // end for j1
		Float_t phys = 0;        
		aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
	    } // end if of threshold condition
	} // for c
    }// end do on pixels
}
//______________________________________________________________________
void AliITSsimulationSPD::SetFluctuations(AliITSpList *pList,Int_t module) {
    //  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
    Float_t  thr=0.0,sigm=0.0;
    Double_t signal,sigma;
    Int_t iz,ix;

    GetThresholds(thr,sigm);
    sigma = (Double_t) sigm;
    for(iz=1;iz<=GetNPixelsZ();iz++){
	for(ix=1;ix<=GetNPixelsX();ix++){
	    signal = sigma*gRandom->Gaus(); 
	    fMapA2->SetHit(iz,ix,signal);
	    // insert in the label-signal-hit list the pixels fired
	    // only by noise
	    pList->AddNoise(iz,ix,module,signal);
	} // end of loop on pixels
    } // end of loop on pixels
}
//______________________________________________________________________
void AliITSsimulationSPD::SetMask() {
    //  Apply a mask to the SPD module. 1% of the pixel channels are
    //  masked. When the database will be ready, the masked pixels
    //  should be read from it.
    Double_t signal;
    Int_t iz,ix,im;
    Float_t totMask;
    Float_t perc = ((AliITSresponseSPD*)fResponse)->GetFractionDead();
    // in this way we get the same set of random numbers for all runs.
    // This is a cluge for now.
    static TRandom *rnd = new TRandom();

    totMask= perc*GetNPixelsZ()*GetNPixelsX();
    for(im=1;im<totMask;im++){
	do{
	    ix=(Int_t)(rnd->Rndm()*(GetNPixelsX()-1.) + 1.);
	} while(ix<=0 || ix>GetNPixelsX());
	do{
	    iz=(Int_t)(rnd->Rndm()*(GetNPixelsZ()-1.) + 1.);
	} while(iz<=0 || iz>GetNPixelsZ());
	signal = -1.;
	fMapA2->SetHit(iz,ix,signal);
    } // end loop on masked pixels
}
//______________________________________________________________________
void AliITSsimulationSPD::CreateHistograms() {
    // Create Histograms
    Int_t i;

    fHis=new TObjArray(GetNPixelsZ());
    for(i=0;i<GetNPixelsZ();i++) {
	TString spdname("spd_");
	Char_t candnum[4];
	sprintf(candnum,"%d",i+1);
	spdname.Append(candnum);
	(*fHis)[i] = new TH1F(spdname.Data(),"SPD maps",
			      GetNPixelsX(),0.,(Float_t) GetNPixelsX());
    } // end for i
}
//______________________________________________________________________
void AliITSsimulationSPD::ResetHistograms() {
    // Reset histograms for this detector
    Int_t i;

    for(i=0;i<GetNPixelsZ();i++ ) {
	if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
    } // end for i
}
//______________________________________________________________________
void AliITSsimulationSPD::WriteSDigits(AliITSpList *pList){
    // Fills the Summable digits Tree
    Int_t i,ni,j,nj;
    static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");

    pList->GetMaxMapIndex(ni,nj);
    for(i=0;i<ni;i++)for(j=0;j<nj;j++){
	if(pList->GetSignalOnly(i,j)>0.0){
	    aliITS->AddSumDigit(*(pList->GetpListItem(i,j)));
//	    cout << "pListSPD: " << *(pList->GetpListItem(i,j)) << endl;
	} // end if
    } // end for i,j
    return;
}
//______________________________________________________________________
void AliITSsimulationSPD::FillMapFrompList(AliITSpList *pList){
    // Fills fMap2A from the pList of Summable digits
    Int_t k,ix;

    for(k=0;k<GetNPixelsZ();k++)for(ix=0;ix<GetNPixelsX();ix++) 
        fMapA2->AddSignal(k,ix,pList->GetSignal(k,ix));
    return;
}