[u/mrichter/AliRoot.git] / TPC / fastSimul / AliTPCclusterFast.cxx

/*
  gSystem->Load("libSTAT.so");
  .x ~/NimStyle.C
 
 .L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
 //
 AliTPCclusterFast::fPRF = new TF1("fprf","gausn");
 AliTPCclusterFast::fTRF = new TF1("ftrf","gausn");
 AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
 AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
 //

 AliTPCclusterFast::Simul("aaa.root",5000); 
 gSystem->Load("libSTAT.so");

 TFile f("aaa.root");
 TTree * tree = (TTree*)f.Get("simul");

*/

#include "TObject.h"
#include "TF1.h"
#include "TMath.h"
#include "TRandom.h"
#include "TVectorD.h"
#include "TMatrixD.h"
#include "TTreeStream.h"

class AliTPCclusterFast: public TObject {
public:
  AliTPCclusterFast();
  virtual ~AliTPCclusterFast();
  void SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz);
  void GenerElectrons();
  void Digitize();
  Double_t GetQtot(Float_t thr, Float_t noise);
  Double_t GetQmax(Float_t thr, Float_t noise);
  Double_t GetNsec();
  static void Simul(const char* simul, Int_t npoints);
public:
  Float_t fMNprim;     // mean number of primary electrons
  //                   //electrons part input
  Int_t   fNprim;      // mean number of primary electrons
  Int_t   fNtot;       // total number of primary electrons 
  Float_t fDiff;       // diffusion sigma
  Float_t fY;          // y position 
  Float_t fZ;          // z postion 
  Float_t fAngleY;     // y angle - tan(y)
  Float_t fAngleZ;     // z angle - tan z
  //
  //
  //                   // electron part simul
  TVectorD fSec;       // number of secondary electrons
  TVectorD fPosY;      //! position y for each electron
  TVectorD fPosZ;      //! position z for each electron
  TVectorD fGain;      //! gg for each electron
  //
  TVectorD fStatY;     // stat Y  
  TVectorD fStatZ;     // stat Y
  //
  // digitization part
  //
  TMatrixD fDigits;    // response matrix
  static TF1* fPRF;    // Pad response
  static TF1* fTRF;    // Time response function 
  ClassDef(AliTPCclusterFast,1)  // container for
};
ClassImp(AliTPCclusterFast)

TF1 *AliTPCclusterFast::fPRF=0;
TF1 *AliTPCclusterFast::fTRF=0;


AliTPCclusterFast::AliTPCclusterFast(){
  //
  //
  fDigits.ResizeTo(5,7);
}

AliTPCclusterFast::~AliTPCclusterFast(){
}


void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz){
  //
  //
  fMNprim = mnprim; fDiff = diff;
  fY=y; fZ=z; 
  fAngleY=ky; fAngleZ=kz;
}
Double_t AliTPCclusterFast::GetNsec(){
  //
  // Generate number of secondary electrons
  // copy of procedure implemented in geant
  //
  const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
  const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
  const Double_t W=20.77E-9;
  Float_t RAN = gRandom->Rndm();
  return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
}

void AliTPCclusterFast::GenerElectrons(){
  //
  //
  //
  //
  const Int_t knMax=1000;
  if (fPosY.GetNrows()<knMax){
    fPosY.ResizeTo(knMax);
    fPosZ.ResizeTo(knMax);
    fGain.ResizeTo(knMax);
    fSec.ResizeTo(knMax);
    fStatY.ResizeTo(3);
    fStatZ.ResizeTo(3);
  }
  fNprim = gRandom->Poisson(fMNprim);  //number of primary electrons
  fNtot=0;
  //
  Double_t sumQ=0;
  Double_t sumYQ=0;
  Double_t sumZQ=0;
  Double_t sumY2Q=0;
  Double_t sumZ2Q=0;
  for (Int_t i=0;i<knMax;i++){ 
    fSec[i]=0;
  }
  for (Int_t iprim=0; iprim<fNprim;iprim++){
    Float_t dN   =  GetNsec();
    fSec[iprim]=dN;
    Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
    Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
    for (Int_t isec=0;isec<=dN;isec++){
      //
      //
      Double_t y = gRandom->Gaus(0,fDiff)+yc;
      Double_t z = gRandom->Gaus(0,fDiff)+zc;
      Double_t gg = -TMath::Log(gRandom->Rndm());
      fPosY[fNtot]=y;
      fPosZ[fNtot]=z;
      fGain[fNtot]=gg;
      fNtot++;
      sumQ+=gg;
      sumYQ+=gg*y;
      sumY2Q+=gg*y*y;
      sumZQ+=gg*z;
      sumZ2Q+=gg*z*z;
      if (fNtot>=knMax) break;
    }
    if (fNtot>=knMax) break;
  }
  if (sumQ>0){
    fStatY[0]=sumQ;
    fStatY[1]=sumYQ/sumQ;
    fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
    fStatZ[0]=sumQ;
    fStatZ[1]=sumZQ/sumQ;
    fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
  }
}

void AliTPCclusterFast::Digitize(){
  //
  //
  //
  // 1. Clear digits
  for (Int_t i=0; i<5;i++)
    for (Int_t j=0; j<7;j++){
      fDigits(i,j)=0;
    }
  //
  // Fill digits
  for (Int_t iel = 0; iel<fNtot; iel++){
    for (Int_t di=-2; di<=2;di++)
      for (Int_t dj=-2; dj<=2;dj++){
	Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
	fac*=fGain[iel];
	fDigits(2+di,3+dj)+=fac;
      }
  }
}


void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
  AliTPCclusterFast fast;
  TTreeSRedirector cstream(fname);
  for (Int_t icl=0; icl<npoints; icl++){
    Float_t nprim=(10+40*gRandom->Rndm());
    Float_t diff =0.01 +0.35*gRandom->Rndm();
    Float_t posY = gRandom->Rndm()-0.5;
    Float_t posZ = gRandom->Rndm()-0.5;
    Float_t ky   = 2.0*(gRandom->Rndm()-0.5);
    Float_t kz   = 2.0*(gRandom->Rndm()-0.5);
    fast.SetParam(nprim,diff,posY,posZ,ky,kz);
    fast.GenerElectrons();
    fast.Digitize();
    cstream<<"simul"<<
      "s.="<<&fast<<
      "\n";
  }
}


Double_t AliTPCclusterFast::GetQtot(Float_t thr, Float_t noise){
  //
  //
  //
  Float_t sum =0;
  for (Int_t i=0;i<5;i++)
    for (Int_t j=0;j<5;j++){
      Float_t amp = fDigits(i,j)+gRandom->Gaus()*noise;
      if (amp>thr) sum+=amp;
    } 
  return sum;
}

Double_t AliTPCclusterFast::GetQmax(Float_t thr, Float_t noise){
  //
  //
  //
  Float_t max =0;
  for (Int_t i=0;i<5;i++)
    for (Int_t j=0;j<5;j++){
      Float_t amp = fDigits(i,j)+gRandom->Gaus()*noise;
      if (amp>thr &&amp>max) max=amp;
    } 
  return max;
}


/*
TH2F *hisL = new TH2F("hisL","hisL",10,10,50,100,0,10);

tree->SetAlias("qmaxNormY","sqrt(sqrt(fDiff^2+fAngleY^2/12+0.5^2))/sqrt(0.5)");
tree->SetAlias("qmaxNormZ","sqrt(sqrt(fDiff^2+fAngleZ^2/12+0.5^2))/sqrt(0.5)");
tree->SetAlias("qmaxNorm","qmaxNormY*qmaxNormZ");


*/
Commit	Line	Data
0716cb63	1	/*
0716cb63	2	gSystem->Load("libSTAT.so");
37045aeb	3	.x ~/NimStyle.C
0716cb63	4
0716cb63	5	.L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
37045aeb	6	//
	7	AliTPCclusterFast::fPRF = new TF1("fprf","gausn");
	8	AliTPCclusterFast::fTRF = new TF1("ftrf","gausn");
	9	AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
	10	AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
	11	//
	12
a314e83c	13	AliTPCclusterFast::Simul("aaa.root",5000);
37045aeb	14	gSystem->Load("libSTAT.so");
	15
	16	TFile f("aaa.root");
	17	TTree * tree = (TTree*)f.Get("simul");
0716cb63	18
	19	*/
	20
	21	#include "TObject.h"
37045aeb	22	#include "TF1.h"
0716cb63	23	#include "TMath.h"
	24	#include "TRandom.h"
	25	#include "TVectorD.h"
37045aeb	26	#include "TMatrixD.h"
0716cb63	27	#include "TTreeStream.h"
	28
	29	class AliTPCclusterFast: public TObject {
	30	public:
	31	AliTPCclusterFast();
	32	virtual ~AliTPCclusterFast();
	33	void SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz);
	34	void GenerElectrons();
37045aeb	35	void Digitize();
a314e83c	36	Double_t GetQtot(Float_t thr, Float_t noise);
a314e83c	37	Double_t GetQmax(Float_t thr, Float_t noise);
37045aeb	38	Double_t GetNsec();
0716cb63	39	static void Simul(const char* simul, Int_t npoints);
	40	public:
	41	Float_t fMNprim; // mean number of primary electrons
	42	// //electrons part input
37045aeb	43	Int_t fNprim; // mean number of primary electrons
0716cb63	44	Int_t fNtot; // total number of primary electrons
	45	Float_t fDiff; // diffusion sigma
	46	Float_t fY; // y position
	47	Float_t fZ; // z postion
	48	Float_t fAngleY; // y angle - tan(y)
	49	Float_t fAngleZ; // z angle - tan z
37045aeb	50	//
37045aeb	51	//
0716cb63	52	// // electron part simul
37045aeb	53	TVectorD fSec; // number of secondary electrons
0716cb63	54	TVectorD fPosY; //! position y for each electron
	55	TVectorD fPosZ; //! position z for each electron
	56	TVectorD fGain; //! gg for each electron
37045aeb	57	//
37045aeb	58	TVectorD fStatY; // stat Y
0716cb63	59	TVectorD fStatZ; // stat Y
0716cb63	60	//
37045aeb	61	// digitization part
	62	//
	63	TMatrixD fDigits; // response matrix
	64	static TF1* fPRF; // Pad response
	65	static TF1* fTRF; // Time response function
0716cb63	66	ClassDef(AliTPCclusterFast,1) // container for
	67	};
	68	ClassImp(AliTPCclusterFast)
	69
37045aeb	70	TF1 *AliTPCclusterFast::fPRF=0;
37045aeb	71	TF1 *AliTPCclusterFast::fTRF=0;
0716cb63	72
	73
	74	AliTPCclusterFast::AliTPCclusterFast(){
37045aeb	75	//
	76	//
	77	fDigits.ResizeTo(5,7);
0716cb63	78	}
37045aeb	79
0716cb63	80	AliTPCclusterFast::~AliTPCclusterFast(){
	81	}
	82
	83
	84	void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz){
	85	//
	86	//
	87	fMNprim = mnprim; fDiff = diff;
	88	fY=y; fZ=z;
	89	fAngleY=ky; fAngleZ=kz;
	90	}
37045aeb	91	Double_t AliTPCclusterFast::GetNsec(){
0716cb63	92	//
37045aeb	93	// Generate number of secondary electrons
37045aeb	94	// copy of procedure implemented in geant
0716cb63	95	//
0716cb63	96	const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
	97	const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
	98	const Double_t W=20.77E-9;
37045aeb	99	Float_t RAN = gRandom->Rndm();
	100	return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)(1-RAN)+TMath::Power(EEND,XEXPO)RAN),YEXPO)/W);
	101	}
	102
	103	void AliTPCclusterFast::GenerElectrons(){
0716cb63	104	//
37045aeb	105	//
	106	//
	107	//
	108	const Int_t knMax=1000;
0716cb63	109	if (fPosY.GetNrows()<knMax){
	110	fPosY.ResizeTo(knMax);
	111	fPosZ.ResizeTo(knMax);
	112	fGain.ResizeTo(knMax);
37045aeb	113	fSec.ResizeTo(knMax);
0716cb63	114	fStatY.ResizeTo(3);
	115	fStatZ.ResizeTo(3);
	116	}
37045aeb	117	fNprim = gRandom->Poisson(fMNprim); //number of primary electrons
0716cb63	118	fNtot=0;
	119	//
	120	Double_t sumQ=0;
	121	Double_t sumYQ=0;
	122	Double_t sumZQ=0;
	123	Double_t sumY2Q=0;
	124	Double_t sumZ2Q=0;
37045aeb	125	for (Int_t i=0;i<knMax;i++){
	126	fSec[i]=0;
	127	}
0716cb63	128	for (Int_t iprim=0; iprim<fNprim;iprim++){
37045aeb	129	Float_t dN = GetNsec();
37045aeb	130	fSec[iprim]=dN;
0716cb63	131	Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
0716cb63	132	Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
37045aeb	133	for (Int_t isec=0;isec<=dN;isec++){
0716cb63	134	//
	135	//
	136	Double_t y = gRandom->Gaus(0,fDiff)+yc;
	137	Double_t z = gRandom->Gaus(0,fDiff)+zc;
	138	Double_t gg = -TMath::Log(gRandom->Rndm());
	139	fPosY[fNtot]=y;
	140	fPosZ[fNtot]=z;
	141	fGain[fNtot]=gg;
	142	fNtot++;
	143	sumQ+=gg;
	144	sumYQ+=gg*y;
	145	sumY2Q+=ggyy;
	146	sumZQ+=gg*z;
	147	sumZ2Q+=ggzz;
	148	if (fNtot>=knMax) break;
	149	}
	150	if (fNtot>=knMax) break;
	151	}
	152	if (sumQ>0){
	153	fStatY[0]=sumQ;
	154	fStatY[1]=sumYQ/sumQ;
	155	fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
	156	fStatZ[0]=sumQ;
	157	fStatZ[1]=sumZQ/sumQ;
	158	fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
	159	}
	160	}
	161
37045aeb	162	void AliTPCclusterFast::Digitize(){
	163	//
	164	//
	165	//
	166	// 1. Clear digits
	167	for (Int_t i=0; i<5;i++)
	168	for (Int_t j=0; j<7;j++){
	169	fDigits(i,j)=0;
	170	}
	171	//
	172	// Fill digits
	173	for (Int_t iel = 0; iel<fNtot; iel++){
	174	for (Int_t di=-2; di<=2;di++)
	175	for (Int_t dj=-2; dj<=2;dj++){
	176	Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
	177	fac*=fGain[iel];
	178	fDigits(2+di,3+dj)+=fac;
	179	}
	180	}
	181	}
	182
	183
0716cb63	184
	185	void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
	186	AliTPCclusterFast fast;
	187	TTreeSRedirector cstream(fname);
	188	for (Int_t icl=0; icl<npoints; icl++){
37045aeb	189	Float_t nprim=(10+40*gRandom->Rndm());
a314e83c	190	Float_t diff =0.01 +0.35*gRandom->Rndm();
0716cb63	191	Float_t posY = gRandom->Rndm()-0.5;
0716cb63	192	Float_t posZ = gRandom->Rndm()-0.5;
a314e83c	193	Float_t ky = 2.0*(gRandom->Rndm()-0.5);
a314e83c	194	Float_t kz = 2.0*(gRandom->Rndm()-0.5);
0716cb63	195	fast.SetParam(nprim,diff,posY,posZ,ky,kz);
0716cb63	196	fast.GenerElectrons();
37045aeb	197	fast.Digitize();
0716cb63	198	cstream<<"simul"<<
	199	"s.="<<&fast<<
	200	"\n";
	201	}
	202	}
	203
	204
a314e83c	205	Double_t AliTPCclusterFast::GetQtot(Float_t thr, Float_t noise){
	206	//
	207	//
	208	//
	209	Float_t sum =0;
	210	for (Int_t i=0;i<5;i++)
	211	for (Int_t j=0;j<5;j++){
	212	Float_t amp = fDigits(i,j)+gRandom->Gaus()*noise;
	213	if (amp>thr) sum+=amp;
	214	}
	215	return sum;
	216	}
	217
	218	Double_t AliTPCclusterFast::GetQmax(Float_t thr, Float_t noise){
	219	//
	220	//
	221	//
	222	Float_t max =0;
	223	for (Int_t i=0;i<5;i++)
	224	for (Int_t j=0;j<5;j++){
	225	Float_t amp = fDigits(i,j)+gRandom->Gaus()*noise;
	226	if (amp>thr &&amp>max) max=amp;
	227	}
	228	return max;
	229	}
	230
	231
0716cb63	232	/*
a314e83c	233	TH2F *hisL = new TH2F("hisL","hisL",10,10,50,100,0,10);
	234
	235	tree->SetAlias("qmaxNormY","sqrt(sqrt(fDiff^2+fAngleY^2/12+0.5^2))/sqrt(0.5)");
	236	tree->SetAlias("qmaxNormZ","sqrt(sqrt(fDiff^2+fAngleZ^2/12+0.5^2))/sqrt(0.5)");
	237	tree->SetAlias("qmaxNorm","qmaxNormY*qmaxNormZ");
	238
	239
0716cb63	240
0716cb63	241	*/