[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 "AliLoader.h"
#include "AliRunLoader.h"

#include <stdlib.h>
#include <Riostream.h>
#include <Riostream.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)
{


  AliRunLoader *inRL, *outRL;//in and out Run Loaders
  AliLoader *ingime, *outgime;// in and out ITSLoaders

  outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
  outgime = outRL->GetLoader("STARTLoader");

#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;
  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 ();

    inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
    ingime = inRL->GetLoader("STARTLoader");
    ingime->LoadHits("READ");//probably it is necessary to load them before
    TClonesArray *STARThitsPhotons = START->Photons ();


   TTree *th = ingime->TreeH();
    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:
    outgime->WriteDigits("OVERWRITE");
  }
}


//------------------------------------------------------------------------
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"
ede9aff7	37	#include "AliPDG.h"
88cb7938	38	#include "AliLoader.h"
88cb7938	39	#include "AliRunLoader.h"
ede9aff7	40
ede9aff7	41	#include <stdlib.h>
ef0750c2	42	#include <Riostream.h>
ef0750c2	43	#include <Riostream.h>
ede9aff7	44
	45	ClassImp(AliSTARTDigitizer)
	46
	47	//___________________________________________
	48	AliSTARTDigitizer::AliSTARTDigitizer() :AliDigitizer()
	49	{
	50	// Default ctor - don't use it
	51	;
	52	}
	53
	54	//___________________________________________
	55	AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager)
	56	:AliDigitizer(manager)
	57	{
	58	cout<<"AliSTARTDigitizer::AliSTARTDigitizer"<<endl;
	59	// ctor which should be used
	60	// fDebug =0;
	61	// if (GetDebug()>2)
	62	// cerr<<"AliSTARTDigitizer::AliSTARTDigitizer"
	63	// <<"(AliRunDigitizer* manager) was processed"<<endl;
	64
	65	}
	66
	67	//------------------------------------------------------------------------
	68	AliSTARTDigitizer::~AliSTARTDigitizer()
	69	{
	70	// Destructor
	71
	72
	73	}
	74
	75	//------------------------------------------------------------------------
	76	Bool_t AliSTARTDigitizer::Init()
	77	{
	78	// Initialization
	79	cout<<"AliSTARTDigitizer::Init"<<endl;
	80	return kTRUE;
	81	}
	82
	83
	84	//---------------------------------------------------------------------
	85
	86	void AliSTARTDigitizer::Exec(Option_t* option)
	87	{
	88
88cb7938	89
	90	AliRunLoader inRL, outRL;//in and out Run Loaders
	91	AliLoader ingime, outgime;// in and out ITSLoaders
	92
	93	outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
	94	outgime = outRL->GetLoader("STARTLoader");
	95
ede9aff7	96	#ifdef DEBUG
	97	cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
	98	#endif
	99	//
	100	// From hits to digits
	101	//
e73d68f2	102	Int_t hit, nhits;
	103	Int_t CountEr[13],CountEl[13]; //!!!
	104	Int_t volume,pmt,tr,tl,sumRight;
e73d68f2	105	Float_t timediff,timeav;
ede9aff7	106	Float_t besttimeright,besttimeleft,meanTime;
e73d68f2	107	Int_t bestRightADC,bestLeftADC;
	108	Float_t besttimeleftGaus, besttimerightGaus;
	109	Float_t timeright[13]={13*0};
	110	Float_t timeleft[13]={13*0};
	111	Float_t channelWidth=2.5; //ps
	112	Int_t channelWidthADC=1; //ps
1ab4b835	113	// Int_t thresholdAmpl=10;
e73d68f2	114
	115	ftimeRightTDC = new TArrayI(12);
	116	ftimeLeftTDC = new TArrayI(12);
	117	fRightADC = new TArrayI(12);
	118	fLeftADC = new TArrayI(12);
	119
	120
ede9aff7	121	fHits = new TClonesArray ("AliSTARThit", 1000);
e73d68f2	122	fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000); //!!!
ede9aff7	123	AliSTART START = (AliSTART) gAlice->GetDetector("START");
	124
	125	AliSTARThit *startHit;
e73d68f2	126	AliSTARThitPhoton *startHitPhoton; //!!!
ede9aff7	127	TBranch *brHits=0;
e73d68f2	128	TBranch *brHitPhoton=0;
ede9aff7	129	fdigits= new AliSTARTdigit();
	130
	131	Int_t nFiles=fManager->GetNinputs();
	132	for (Int_t inputFile=0; inputFile<nFiles; inputFile++) {
e73d68f2	133
ede9aff7	134
e73d68f2	135
ede9aff7	136	besttimeright=9999.;
	137	besttimeleft=9999.;
	138	Int_t timeDiff=0;
	139	Int_t timeAv=0;
e73d68f2	140	sumRight=0;
	141	for (Int_t i0=0; i0<13; i0++)
	142	{
	143	timeright[i0]=0; timeleft[i0]=0;
	144	CountEr[i0]=0; CountEl[i0]=0;
	145	}
ede9aff7	146	TClonesArray *STARThits = START->Hits ();
88cb7938	147
	148	inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
	149	ingime = inRL->GetLoader("STARTLoader");
	150	ingime->LoadHits("READ");//probably it is necessary to load them before
e73d68f2	151	TClonesArray *STARThitsPhotons = START->Photons ();
ede9aff7	152
88cb7938	153
88cb7938	154	TTree *th = ingime->TreeH();
ede9aff7	155	brHits = th->GetBranch("START");
e73d68f2	156	brHitPhoton = th->GetBranch("STARThitPhoton");
ede9aff7	157	if (brHits) {
e73d68f2	158	START->SetHitsAddressBranch(brHits,brHitPhoton);
ede9aff7	159	}else{
00907af4	160	cerr<<"EXEC Branch START hit not found"<<endl;
00907af4	161	exit(111);
ede9aff7	162	}
ede9aff7	163	Int_t ntracks = (Int_t) th->GetEntries();
e73d68f2	164	cout<<" ntracks "<<ntracks<<endl;
ede9aff7	165	if (ntracks<=0) return;
e73d68f2	166	// Start loop on tracks in the photon hits containers
	167	// for amplitude
	168	if(brHitPhoton) {
	169	for (Int_t track=0; track<ntracks;track++) {
	170	brHitPhoton -> GetEntry(track);;
	171	nhits = STARThitsPhotons->GetEntriesFast();
	172	for (hit=0;hit<nhits;hit++) {
	173	startHitPhoton = (AliSTARThitPhoton*)
	174	STARThitsPhotons ->UncheckedAt(hit);
	175	pmt=startHitPhoton->fPmt;
	176	volume = startHitPhoton->fArray;
	177	if(RegisterPhotoE(startHitPhoton))
	178	{
	179	if (volume == 1) CountEr[pmt]++;
	180	if (volume == 2) CountEl[pmt]++;
	181	}
	182	} //hit photons
	183	} //track photons
	184	} // was photons
	185
ede9aff7	186	// Start loop on tracks in the hits containers
	187	for (Int_t track=0; track<ntracks;track++) {
	188	brHits->GetEntry(track);
	189	nhits = STARThits->GetEntriesFast();
	190	for (hit=0;hit<nhits;hit++) {
	191	startHit = (AliSTARThit*) STARThits->UncheckedAt(hit);
	192	pmt=startHit->fPmt;
	193	volume = startHit->fVolume;
	194	if(volume==1){
e73d68f2	195	timeright[pmt] = startHit->fTime;
	196	if(timeright[pmt]<besttimeright)
	197	//&&CountEr[pmt-1]>thresholdAmpl)
	198	{
	199	besttimeright=timeright[pmt];
ede9aff7	200	} //timeright
	201	}//time for right shoulder
	202	if(volume==2){
e73d68f2	203	timeleft[pmt] = startHit->fTime;
	204	if(timeleft[pmt]<besttimeleft)
	205	//&&CountEl[pmt-1]>thresholdAmpl)
	206	{
	207	besttimeleft=timeleft[pmt];
	208
ede9aff7	209	} //timeleftbest
	210	}//time for left shoulder
	211	} //hit loop
	212	} //track loop
	213
e73d68f2	214	// z position
	215	cout<<" right time "<<besttimeright<<
	216	" right distance "<<besttimeright*30<<endl;;
	217	cout<<" left time "<<besttimeleft<<
	218	" right distance "<<besttimeleft*30<<endl;;
	219
	220
ede9aff7	221	//folding with experimental time distribution
e73d68f2	222
	223	besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
	224	// cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
1ab4b835	225	bestRightADC=Int_t (besttimerightGaus*1000/channelWidth);
7115262b	226	Float_t koef=69.7/350.;
7115262b	227	besttimeleft=koef*besttimeleft;
e73d68f2	228	besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
e73d68f2	229
1ab4b835	230	bestLeftADC=Int_t (besttimeleftGaus*1000/channelWidth);
ede9aff7	231	timediff=besttimerightGaus-besttimeleftGaus;
e73d68f2	232	cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
ede9aff7	233	meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
e73d68f2	234	if ( TMath::Abs(timediff)<TMath::Abs(0.3) )
	235	{
	236	Float_t t1=1000.*besttimeleftGaus;
	237	Float_t t2=1000.*besttimerightGaus;
	238	t1=t1/channelWidth; //time in ps to channelWidth
	239	t2=t2/channelWidth; //time in ps to channelWidth
	240	timeav=(t1+t2)/2.;
	241
	242	// Time to TDC signal
	243	// 256 channels for timediff, range 1ns
	244
	245	timediff=512+1000*timediff/channelWidth; // time in ps
	246
	247	timeAv = (Int_t)(timeav); // time channel numbres
	248	timeDiff = (Int_t)(timediff); // time channel numbres
	249	// fill digits
	250	fdigits->SetTimeBestLeft(bestLeftADC);
	251	fdigits->SetTimeBestRight(bestRightADC);
	252	fdigits->SetMeanTime(timeAv);
	253	fdigits->SetTimeDiff(timeDiff);
	254	for (Int_t i=0; i<12; i++)
	255	{
	256	// fill TDC
	257	timeright[i+1]=gRandom->Gaus(timeright[i+1],0.05);
	258	timeleft[i+1]=gRandom->Gaus(timeleft[i+1],0.05);
	259	tr= Int_t (timeright[i+1]*1000/channelWidth);
	260	if(tr<200) tr=0;
	261	tl= Int_t (timeleft[i+1]*1000/channelWidth);
	262	if(tl<1000) tl=0;
	263
	264	ftimeRightTDC->AddAt(tr,i);
	265	ftimeLeftTDC->AddAt(tl,i);
	266	//fill ADC
	267	Int_t al=( Int_t ) CountEl[i+1]/ channelWidthADC;
	268	Int_t ar=( Int_t ) CountEr[i+1]/ channelWidthADC;
	269	fRightADC->AddAt(ar,i);
	270	fLeftADC ->AddAt(al,i);
	271	sumRight+=CountEr[i+1];
	272	}
	273	fdigits->SetTimeRight(*ftimeRightTDC);
	274	fdigits->SetTimeLeft(*ftimeLeftTDC);
	275	fdigits->SetADCRight(*fRightADC);
	276	fdigits->SetADCLeft(*fLeftADC);
	277	// cout<<" before sum"<<endl;
	278	fdigits->SetSumADCRight(sumRight);
	279	}
ede9aff7	280	else
ede9aff7	281	{timeAv=999999; timeDiff=99999;}
88cb7938	282
	283	// trick to find out output dir:
	284	outgime->WriteDigits("OVERWRITE");
ede9aff7	285	}
ede9aff7	286	}
	287
	288
e73d68f2	289	//------------------------------------------------------------------------
	290	Bool_t AliSTARTDigitizer::RegisterPhotoE(AliSTARThitPhoton *hit)
	291	{
	292	Double_t P = 0.2;
	293	Double_t p;
	294
	295	p = gRandom->Rndm();
	296	if (p > P)
	297	return kFALSE;
	298
	299	return kTRUE;
	300	}