new Digits structure
[u/mrichter/AliRoot.git] / START / AliSTARTDigitizer.cxx
CommitLineData
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>
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"
33#include "AliRunDigitizer.h"
e73d68f2 34#include "AliHeader.h"
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
43ClassImp(AliSTARTDigitizer)
44
45//___________________________________________
46 AliSTARTDigitizer::AliSTARTDigitizer() :AliDigitizer()
47{
48// Default ctor - don't use it
49 ;
50}
51
52//___________________________________________
53AliSTARTDigitizer::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//------------------------------------------------------------------------
66AliSTARTDigitizer::~AliSTARTDigitizer()
67{
68// Destructor
69
70
71}
72
73 //------------------------------------------------------------------------
74Bool_t AliSTARTDigitizer::Init()
75{
76// Initialization
77 cout<<"AliSTARTDigitizer::Init"<<endl;
78 return kTRUE;
79}
80
81
82//---------------------------------------------------------------------
83
84void 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;
149 exit(111);
ede9aff7 150 }
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.;
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
269 {timeAv=999999; timeDiff=99999;}
e73d68f2 270
271 // trick to find out output dir:
ede9aff7 272 TTree *outTree = fManager->GetTreeD();
273 if (!outTree) {
00907af4 274 cerr<<"something wrong with output...."<<endl;
275 exit(111);
ede9aff7 276 }
e73d68f2 277
ede9aff7 278 TDirectory *wd = gDirectory;
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();
284 }
ede9aff7 285}
286
287
e73d68f2 288//------------------------------------------------------------------------
289Bool_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}