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