Transition to NewIO
[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"
88cb7938 38#include "AliLoader.h"
39#include "AliRunLoader.h"
ede9aff7 40
41#include <stdlib.h>
ef0750c2 42#include <Riostream.h>
43#include <Riostream.h>
ede9aff7 44
45ClassImp(AliSTARTDigitizer)
46
47//___________________________________________
48 AliSTARTDigitizer::AliSTARTDigitizer() :AliDigitizer()
49{
50// Default ctor - don't use it
51 ;
52}
53
54//___________________________________________
55AliSTARTDigitizer::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//------------------------------------------------------------------------
68AliSTARTDigitizer::~AliSTARTDigitizer()
69{
70// Destructor
71
72
73}
74
75 //------------------------------------------------------------------------
76Bool_t AliSTARTDigitizer::Init()
77{
78// Initialization
79 cout<<"AliSTARTDigitizer::Init"<<endl;
80 return kTRUE;
81}
82
83
84//---------------------------------------------------------------------
85
86void 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
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;
161 exit(111);
ede9aff7 162 }
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.;
227 besttimeleft=koef*besttimeleft;
e73d68f2 228 besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
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
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//------------------------------------------------------------------------
290Bool_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}