[u/mrichter/AliRoot.git] / START / AliSTARTDigitizer.cxx


/**************************************************************************
 * Copyright(c) 1998-2000, 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.                  *
 **************************************************************************/


#include <TTree.h> 
#include <TVector.h>
#include <TObjArray.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TRandom.h>
#include <TArrayI.h>
#include <TH1.h>


#include "AliSTARTDigitizer.h"
#include "AliSTART.h"
#include "AliSTARThit.h"
#include "AliSTARThitPhoton.h"
#include "AliSTARTdigit.h"
#include "AliRunDigitizer.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliRun.h"
#include "AliPDG.h"

#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>

ClassImp(AliSTARTDigitizer)

//___________________________________________
  AliSTARTDigitizer::AliSTARTDigitizer()  :AliDigitizer()
{
// Default ctor - don't use it
  ;
}

//___________________________________________
AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager) 
    :AliDigitizer(manager) 
{
	cout<<"AliSTARTDigitizer::AliSTARTDigitizer"<<endl;
// ctor which should be used
//  fDebug =0;
 // if (GetDebug()>2)
  //  cerr<<"AliSTARTDigitizer::AliSTARTDigitizer"
   //     <<"(AliRunDigitizer* manager) was processed"<<endl;

}

//------------------------------------------------------------------------
AliSTARTDigitizer::~AliSTARTDigitizer()
{
// Destructor


}

 //------------------------------------------------------------------------
Bool_t AliSTARTDigitizer::Init()
{
// Initialization
 cout<<"AliSTARTDigitizer::Init"<<endl;
 return kTRUE;
}
 

//---------------------------------------------------------------------

void AliSTARTDigitizer::Exec(Option_t* option)
{

  cout<<" AliSTARTDigitizer::Exec"<<endl;
#ifdef DEBUG
  cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
#endif
  //
  // From hits to digits
  //
  Int_t hit, nhits;
  Int_t CountEr[13],CountEl[13];							//!!!
  Int_t volume,pmt,tr,tl,sumRight;
  char nameDigits[20];
  Float_t timediff,timeav;
  Float_t besttimeright,besttimeleft,meanTime;
  Int_t  bestRightADC,bestLeftADC;
  Float_t besttimeleftGaus, besttimerightGaus;
  Float_t timeright[13]={13*0};
  Float_t timeleft[13]={13*0};
  Float_t channelWidth=2.5; //ps
  Int_t channelWidthADC=1; //ps
  Int_t thresholdAmpl=10;

  ftimeRightTDC = new TArrayI(12); 
  ftimeLeftTDC = new TArrayI(12); 
  fRightADC = new TArrayI(12); 
  fLeftADC = new TArrayI(12); 
  

  fHits = new TClonesArray ("AliSTARThit", 1000);
  fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000);			//!!!
  AliSTART *START  = (AliSTART*) gAlice->GetDetector("START");

  AliSTARThit  *startHit;
  AliSTARThitPhoton  *startHitPhoton;						//!!!
  TBranch *brHits=0;
  TBranch *brHitPhoton=0;
  fdigits= new AliSTARTdigit();

  Int_t nFiles=fManager->GetNinputs();
  for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {

 
    besttimeright=9999.;
    besttimeleft=9999.;
    Int_t timeDiff=0;
    Int_t timeAv=0;
    sumRight=0;
    for (Int_t i0=0; i0<13; i0++)
      {
	timeright[i0]=0; timeleft[i0]=0;
	CountEr[i0]=0;   CountEl[i0]=0;
      }
    TClonesArray *STARThits = START->Hits ();
    TClonesArray *STARThitsPhotons = START->Photons ();

   TTree *th = fManager->GetInputTreeH(inputFile);
    brHits = th->GetBranch("START");
    brHitPhoton = th->GetBranch("STARThitPhoton");
    if (brHits) {
      START->SetHitsAddressBranch(brHits,brHitPhoton);
    }else{
      cerr<<"EXEC Branch START hit not found"<<endl;
      exit(111);
    } 
    Int_t ntracks    = (Int_t) th->GetEntries();
    cout<<" ntracks "<<ntracks<<endl;
    if (ntracks<=0) return;
    // Start loop on tracks in the photon hits containers 
    // for amplitude
    if(brHitPhoton) {
      for (Int_t track=0; track<ntracks;track++) {
	brHitPhoton -> GetEntry(track);;
	nhits = STARThitsPhotons->GetEntriesFast();
	for (hit=0;hit<nhits;hit++) {
	  startHitPhoton   = (AliSTARThitPhoton*) 
	   STARThitsPhotons ->UncheckedAt(hit);
	  pmt=startHitPhoton->fPmt;
	  volume = startHitPhoton->fArray;
	  if(RegisterPhotoE(startHitPhoton)) 
	    {
	      if (volume == 1) CountEr[pmt]++;
	      if (volume == 2) CountEl[pmt]++;
	    }
	} //hit photons
      } //track photons
    } // was photons

    // Start loop on tracks in the hits containers
    for (Int_t track=0; track<ntracks;track++) {
      brHits->GetEntry(track);
      nhits = STARThits->GetEntriesFast();
      for (hit=0;hit<nhits;hit++) {
	startHit   = (AliSTARThit*) STARThits->UncheckedAt(hit);
	pmt=startHit->fPmt;
	volume = startHit->fVolume;
	if(volume==1){
	  timeright[pmt] = startHit->fTime;
	  if(timeright[pmt]<besttimeright)
	    //&&CountEr[pmt-1]>thresholdAmpl)
	    {
	    besttimeright=timeright[pmt];
	  } //timeright
	}//time for right shoulder
	if(volume==2){            
	  timeleft[pmt] = startHit->fTime;
	  if(timeleft[pmt]<besttimeleft)
	    //&&CountEl[pmt-1]>thresholdAmpl) 
	    {
	    besttimeleft=timeleft[pmt];
	    
	  } //timeleftbest
	}//time for left shoulder
      } //hit loop
    } //track loop
 
    // z position
    cout<<" right time  "<<besttimeright<<
      " right distance "<<besttimeright*30<<endl;;
    cout<<" left time  "<<besttimeleft<<
      " right distance "<<besttimeleft*30<<endl;;
  

    //folding with experimental time distribution
    
    besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
    //    cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
    bestRightADC=(Int_t) besttimerightGaus*1000/channelWidth;
    Float_t koef=69.7/350.;
    besttimeleft=koef*besttimeleft;
    besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
    
    bestLeftADC=(Int_t) besttimeleftGaus*1000/channelWidth;
    timediff=besttimerightGaus-besttimeleftGaus;
    cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
    meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
    if ( TMath::Abs(timediff)<TMath::Abs(0.3) ) 
      {
	Float_t t1=1000.*besttimeleftGaus;
	Float_t t2=1000.*besttimerightGaus;
	t1=t1/channelWidth;   //time in ps to channelWidth
	t2=t2/channelWidth;   //time in ps to channelWidth
	timeav=(t1+t2)/2.;
	
	// Time to TDC signal
	// 256 channels for timediff, range 1ns
	
	timediff=512+1000*timediff/channelWidth; // time in ps 
	
	timeAv = (Int_t)(timeav);   // time  channel numbres
	timeDiff = (Int_t)(timediff); // time  channel numbres
	//       fill digits
	fdigits->SetTimeBestLeft(bestLeftADC);
	fdigits->SetTimeBestRight(bestRightADC);
	fdigits->SetMeanTime(timeAv);
	fdigits->SetTimeDiff(timeDiff);
	for (Int_t i=0; i<12; i++)
	  {
	    //  fill TDC
	    timeright[i+1]=gRandom->Gaus(timeright[i+1],0.05);
	    timeleft[i+1]=gRandom->Gaus(timeleft[i+1],0.05);
	    tr= Int_t (timeright[i+1]*1000/channelWidth); 
	    if(tr<200) tr=0;
	    tl= Int_t (timeleft[i+1]*1000/channelWidth); 
	    if(tl<1000) tl=0;
	    
	    ftimeRightTDC->AddAt(tr,i);
	    ftimeLeftTDC->AddAt(tl,i);
	    //fill ADC
	    Int_t al=( Int_t ) CountEl[i+1]/ channelWidthADC;
	    Int_t ar=( Int_t ) CountEr[i+1]/ channelWidthADC;
	    fRightADC->AddAt(ar,i);
	    fLeftADC ->AddAt(al,i);
	    sumRight+=CountEr[i+1];
	  }
	fdigits->SetTimeRight(*ftimeRightTDC);
	fdigits->SetTimeLeft(*ftimeLeftTDC);
	fdigits->SetADCRight(*fRightADC);
	fdigits->SetADCLeft(*fLeftADC);
	// cout<<" before sum"<<endl;
	fdigits->SetSumADCRight(sumRight);
      }
    else
      {timeAv=999999; timeDiff=99999;}
    
    // trick to find out output dir:
    TTree *outTree = fManager->GetTreeD();
    if (!outTree) {
      cerr<<"something wrong with output...."<<endl;
      exit(111);
    }

    TDirectory *wd = gDirectory;
    outTree->GetDirectory()->cd();
    sprintf(nameDigits,"START_D_%d",fManager->GetOutputEventNr());
    fdigits->Write(nameDigits);
    cout<<nameDigits<<endl;
    wd->cd();
  }
}


//------------------------------------------------------------------------
Bool_t AliSTARTDigitizer::RegisterPhotoE(AliSTARThitPhoton *hit)
{
    Double_t    P = 0.2;    
    Double_t    p;
    
    p = gRandom->Rndm();
    if (p > P)
      return kFALSE;
    
    return kTRUE;
}
Commit	Line	Data
ede9aff7	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17
	18	#include <TTree.h>
	19	#include <TVector.h>
	20	#include <TObjArray.h>
	21	#include <TFile.h>
	22	#include <TDirectory.h>
	23	#include <TRandom.h>
e73d68f2	24	#include <TArrayI.h>
e73d68f2	25	#include <TH1.h>
ede9aff7	26
	27
	28	#include "AliSTARTDigitizer.h"
	29	#include "AliSTART.h"
	30	#include "AliSTARThit.h"
e73d68f2	31	#include "AliSTARThitPhoton.h"
ede9aff7	32	#include "AliSTARTdigit.h"
ede9aff7	33	#include "AliRunDigitizer.h"
e73d68f2	34	#include "AliHeader.h"
e73d68f2	35	#include "AliGenEventHeader.h"
ede9aff7	36	#include "AliRun.h"
	37	#include "AliPDG.h"
	38
	39	#include <stdlib.h>
	40	#include <iostream.h>
	41	#include <fstream.h>
	42
	43	ClassImp(AliSTARTDigitizer)
	44
	45	//___________________________________________
	46	AliSTARTDigitizer::AliSTARTDigitizer() :AliDigitizer()
	47	{
	48	// Default ctor - don't use it
	49	;
	50	}
	51
	52	//___________________________________________
	53	AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager)
	54	:AliDigitizer(manager)
	55	{
	56	cout<<"AliSTARTDigitizer::AliSTARTDigitizer"<<endl;
	57	// ctor which should be used
	58	// fDebug =0;
	59	// if (GetDebug()>2)
	60	// cerr<<"AliSTARTDigitizer::AliSTARTDigitizer"
	61	// <<"(AliRunDigitizer* manager) was processed"<<endl;
	62
	63	}
	64
	65	//------------------------------------------------------------------------
	66	AliSTARTDigitizer::~AliSTARTDigitizer()
	67	{
	68	// Destructor
	69
	70
	71	}
	72
	73	//------------------------------------------------------------------------
	74	Bool_t AliSTARTDigitizer::Init()
	75	{
	76	// Initialization
	77	cout<<"AliSTARTDigitizer::Init"<<endl;
	78	return kTRUE;
	79	}
	80
	81
	82	//---------------------------------------------------------------------
	83
	84	void AliSTARTDigitizer::Exec(Option_t* option)
	85	{
	86
e73d68f2	87	cout<<" AliSTARTDigitizer::Exec"<<endl;
ede9aff7	88	#ifdef DEBUG
	89	cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
	90	#endif
	91	//
	92	// From hits to digits
	93	//
e73d68f2	94	Int_t hit, nhits;
	95	Int_t CountEr[13],CountEl[13]; //!!!
	96	Int_t volume,pmt,tr,tl,sumRight;
ede9aff7	97	char nameDigits[20];
e73d68f2	98	Float_t timediff,timeav;
ede9aff7	99	Float_t besttimeright,besttimeleft,meanTime;
e73d68f2	100	Int_t bestRightADC,bestLeftADC;
	101	Float_t besttimeleftGaus, besttimerightGaus;
	102	Float_t timeright[13]={13*0};
	103	Float_t timeleft[13]={13*0};
	104	Float_t channelWidth=2.5; //ps
	105	Int_t channelWidthADC=1; //ps
	106	Int_t thresholdAmpl=10;
	107
	108	ftimeRightTDC = new TArrayI(12);
	109	ftimeLeftTDC = new TArrayI(12);
	110	fRightADC = new TArrayI(12);
	111	fLeftADC = new TArrayI(12);
	112
	113
ede9aff7	114	fHits = new TClonesArray ("AliSTARThit", 1000);
e73d68f2	115	fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000); //!!!
ede9aff7	116	AliSTART START = (AliSTART) gAlice->GetDetector("START");
	117
	118	AliSTARThit *startHit;
e73d68f2	119	AliSTARThitPhoton *startHitPhoton; //!!!
ede9aff7	120	TBranch *brHits=0;
e73d68f2	121	TBranch *brHitPhoton=0;
ede9aff7	122	fdigits= new AliSTARTdigit();
	123
	124	Int_t nFiles=fManager->GetNinputs();
	125	for (Int_t inputFile=0; inputFile<nFiles; inputFile++) {
e73d68f2	126
ede9aff7	127
e73d68f2	128
ede9aff7	129	besttimeright=9999.;
	130	besttimeleft=9999.;
	131	Int_t timeDiff=0;
	132	Int_t timeAv=0;
e73d68f2	133	sumRight=0;
	134	for (Int_t i0=0; i0<13; i0++)
	135	{
	136	timeright[i0]=0; timeleft[i0]=0;
	137	CountEr[i0]=0; CountEl[i0]=0;
	138	}
ede9aff7	139	TClonesArray *STARThits = START->Hits ();
e73d68f2	140	TClonesArray *STARThitsPhotons = START->Photons ();
ede9aff7	141
	142	TTree *th = fManager->GetInputTreeH(inputFile);
	143	brHits = th->GetBranch("START");
e73d68f2	144	brHitPhoton = th->GetBranch("STARThitPhoton");
ede9aff7	145	if (brHits) {
e73d68f2	146	START->SetHitsAddressBranch(brHits,brHitPhoton);
ede9aff7	147	}else{
00907af4	148	cerr<<"EXEC Branch START hit not found"<<endl;
00907af4	149	exit(111);
ede9aff7	150	}
ede9aff7	151	Int_t ntracks = (Int_t) th->GetEntries();
e73d68f2	152	cout<<" ntracks "<<ntracks<<endl;
ede9aff7	153	if (ntracks<=0) return;
e73d68f2	154	// Start loop on tracks in the photon hits containers
	155	// for amplitude
	156	if(brHitPhoton) {
	157	for (Int_t track=0; track<ntracks;track++) {
	158	brHitPhoton -> GetEntry(track);;
	159	nhits = STARThitsPhotons->GetEntriesFast();
	160	for (hit=0;hit<nhits;hit++) {
	161	startHitPhoton = (AliSTARThitPhoton*)
	162	STARThitsPhotons ->UncheckedAt(hit);
	163	pmt=startHitPhoton->fPmt;
	164	volume = startHitPhoton->fArray;
	165	if(RegisterPhotoE(startHitPhoton))
	166	{
	167	if (volume == 1) CountEr[pmt]++;
	168	if (volume == 2) CountEl[pmt]++;
	169	}
	170	} //hit photons
	171	} //track photons
	172	} // was photons
	173
ede9aff7	174	// Start loop on tracks in the hits containers
	175	for (Int_t track=0; track<ntracks;track++) {
	176	brHits->GetEntry(track);
	177	nhits = STARThits->GetEntriesFast();
	178	for (hit=0;hit<nhits;hit++) {
	179	startHit = (AliSTARThit*) STARThits->UncheckedAt(hit);
	180	pmt=startHit->fPmt;
	181	volume = startHit->fVolume;
	182	if(volume==1){
e73d68f2	183	timeright[pmt] = startHit->fTime;
	184	if(timeright[pmt]<besttimeright)
	185	//&&CountEr[pmt-1]>thresholdAmpl)
	186	{
	187	besttimeright=timeright[pmt];
ede9aff7	188	} //timeright
	189	}//time for right shoulder
	190	if(volume==2){
e73d68f2	191	timeleft[pmt] = startHit->fTime;
	192	if(timeleft[pmt]<besttimeleft)
	193	//&&CountEl[pmt-1]>thresholdAmpl)
	194	{
	195	besttimeleft=timeleft[pmt];
	196
ede9aff7	197	} //timeleftbest
	198	}//time for left shoulder
	199	} //hit loop
	200	} //track loop
	201
e73d68f2	202	// z position
	203	cout<<" right time "<<besttimeright<<
	204	" right distance "<<besttimeright*30<<endl;;
	205	cout<<" left time "<<besttimeleft<<
	206	" right distance "<<besttimeleft*30<<endl;;
	207
	208
ede9aff7	209	//folding with experimental time distribution
e73d68f2	210
	211	besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
	212	// cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
	213	bestRightADC=(Int_t) besttimerightGaus*1000/channelWidth;
7115262b	214	Float_t koef=69.7/350.;
7115262b	215	besttimeleft=koef*besttimeleft;
e73d68f2	216	besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
	217
	218	bestLeftADC=(Int_t) besttimeleftGaus*1000/channelWidth;
ede9aff7	219	timediff=besttimerightGaus-besttimeleftGaus;
e73d68f2	220	cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
ede9aff7	221	meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
e73d68f2	222	if ( TMath::Abs(timediff)<TMath::Abs(0.3) )
	223	{
	224	Float_t t1=1000.*besttimeleftGaus;
	225	Float_t t2=1000.*besttimerightGaus;
	226	t1=t1/channelWidth; //time in ps to channelWidth
	227	t2=t2/channelWidth; //time in ps to channelWidth
	228	timeav=(t1+t2)/2.;
	229
	230	// Time to TDC signal
	231	// 256 channels for timediff, range 1ns
	232
	233	timediff=512+1000*timediff/channelWidth; // time in ps
	234
	235	timeAv = (Int_t)(timeav); // time channel numbres
	236	timeDiff = (Int_t)(timediff); // time channel numbres
	237	// fill digits
	238	fdigits->SetTimeBestLeft(bestLeftADC);
	239	fdigits->SetTimeBestRight(bestRightADC);
	240	fdigits->SetMeanTime(timeAv);
	241	fdigits->SetTimeDiff(timeDiff);
	242	for (Int_t i=0; i<12; i++)
	243	{
	244	// fill TDC
	245	timeright[i+1]=gRandom->Gaus(timeright[i+1],0.05);
	246	timeleft[i+1]=gRandom->Gaus(timeleft[i+1],0.05);
	247	tr= Int_t (timeright[i+1]*1000/channelWidth);
	248	if(tr<200) tr=0;
	249	tl= Int_t (timeleft[i+1]*1000/channelWidth);
	250	if(tl<1000) tl=0;
	251
	252	ftimeRightTDC->AddAt(tr,i);
	253	ftimeLeftTDC->AddAt(tl,i);
	254	//fill ADC
	255	Int_t al=( Int_t ) CountEl[i+1]/ channelWidthADC;
	256	Int_t ar=( Int_t ) CountEr[i+1]/ channelWidthADC;
	257	fRightADC->AddAt(ar,i);
	258	fLeftADC ->AddAt(al,i);
	259	sumRight+=CountEr[i+1];
	260	}
	261	fdigits->SetTimeRight(*ftimeRightTDC);
	262	fdigits->SetTimeLeft(*ftimeLeftTDC);
	263	fdigits->SetADCRight(*fRightADC);
	264	fdigits->SetADCLeft(*fLeftADC);
	265	// cout<<" before sum"<<endl;
	266	fdigits->SetSumADCRight(sumRight);
	267	}
ede9aff7	268	else
ede9aff7	269	{timeAv=999999; timeDiff=99999;}
e73d68f2	270
e73d68f2	271	// trick to find out output dir:
ede9aff7	272	TTree *outTree = fManager->GetTreeD();
ede9aff7	273	if (!outTree) {
00907af4	274	cerr<<"something wrong with output...."<<endl;
00907af4	275	exit(111);
ede9aff7	276	}
e73d68f2	277
ede9aff7	278	TDirectory *wd = gDirectory;
ede9aff7	279	outTree->GetDirectory()->cd();
e73d68f2	280	sprintf(nameDigits,"START_D_%d",fManager->GetOutputEventNr());
ede9aff7	281	fdigits->Write(nameDigits);
e73d68f2	282	cout<<nameDigits<<endl;
ede9aff7	283	wd->cd();
ede9aff7	284	}
ede9aff7	285	}
	286
	287
e73d68f2	288	//------------------------------------------------------------------------
	289	Bool_t AliSTARTDigitizer::RegisterPhotoE(AliSTARThitPhoton *hit)
	290	{
	291	Double_t P = 0.2;
	292	Double_t p;
	293
	294	p = gRandom->Rndm();
	295	if (p > P)
	296	return kFALSE;
	297
	298	return kTRUE;
	299	}