90fd8be006652d2fb61161f7f5a57d95dce7ecc3
[u/mrichter/AliRoot.git] / START / AliSTARTDigitizer.cxx
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>
24 #include <TArrayI.h>
25 #include <TH1.h>
26
27
28 #include "AliSTARTDigitizer.h"
29 #include "AliSTART.h"
30 #include "AliSTARThit.h"
31 #include "AliSTARThitPhoton.h"
32 #include "AliSTARTdigit.h"
33 #include "AliRunDigitizer.h"
34 #include "AliHeader.h"
35 #include "AliGenEventHeader.h"
36 #include "AliRun.h"
37 #include "AliPDG.h"
38
39 #include <stdlib.h>
40 #include <Riostream.h>
41 #include <Riostream.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
87   cout<<" AliSTARTDigitizer::Exec"<<endl;
88 #ifdef DEBUG
89   cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
90 #endif
91   //
92   // From hits to digits
93   //
94   Int_t hit, nhits;
95   Int_t CountEr[13],CountEl[13];                                                        //!!!
96   Int_t volume,pmt,tr,tl,sumRight;
97   char nameDigits[20];
98   Float_t timediff,timeav;
99   Float_t besttimeright,besttimeleft,meanTime;
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
114   fHits = new TClonesArray ("AliSTARThit", 1000);
115   fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000);                     //!!!
116   AliSTART *START  = (AliSTART*) gAlice->GetDetector("START");
117
118   AliSTARThit  *startHit;
119   AliSTARThitPhoton  *startHitPhoton;                                           //!!!
120   TBranch *brHits=0;
121   TBranch *brHitPhoton=0;
122   fdigits= new AliSTARTdigit();
123
124   Int_t nFiles=fManager->GetNinputs();
125   for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {
126
127  
128
129     besttimeright=9999.;
130     besttimeleft=9999.;
131     Int_t timeDiff=0;
132     Int_t timeAv=0;
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       }
139     TClonesArray *STARThits = START->Hits ();
140     TClonesArray *STARThitsPhotons = START->Photons ();
141
142    TTree *th = fManager->GetInputTreeH(inputFile);
143     brHits = th->GetBranch("START");
144     brHitPhoton = th->GetBranch("STARThitPhoton");
145     if (brHits) {
146       START->SetHitsAddressBranch(brHits,brHitPhoton);
147     }else{
148       cerr<<"EXEC Branch START hit not found"<<endl;
149       exit(111);
150     } 
151     Int_t ntracks    = (Int_t) th->GetEntries();
152     cout<<" ntracks "<<ntracks<<endl;
153     if (ntracks<=0) return;
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
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){
183           timeright[pmt] = startHit->fTime;
184           if(timeright[pmt]<besttimeright)
185             //&&CountEr[pmt-1]>thresholdAmpl)
186             {
187             besttimeright=timeright[pmt];
188           } //timeright
189         }//time for right shoulder
190         if(volume==2){            
191           timeleft[pmt] = startHit->fTime;
192           if(timeleft[pmt]<besttimeleft)
193             //&&CountEl[pmt-1]>thresholdAmpl) 
194             {
195             besttimeleft=timeleft[pmt];
196             
197           } //timeleftbest
198         }//time for left shoulder
199       } //hit loop
200     } //track loop
201  
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
209     //folding with experimental time distribution
210     
211     besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
212     //    cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
213     bestRightADC=Int_t (besttimerightGaus*1000/channelWidth);
214     Float_t koef=69.7/350.;
215     besttimeleft=koef*besttimeleft;
216     besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
217     
218     bestLeftADC=Int_t (besttimeleftGaus*1000/channelWidth);
219     timediff=besttimerightGaus-besttimeleftGaus;
220     cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
221     meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
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       }
268     else
269       {timeAv=999999; timeDiff=99999;}
270     
271     // trick to find out output dir:
272     TTree *outTree = fManager->GetTreeD();
273     if (!outTree) {
274       cerr<<"something wrong with output...."<<endl;
275       exit(111);
276     }
277
278     TDirectory *wd = gDirectory;
279     outTree->GetDirectory()->cd();
280     sprintf(nameDigits,"START_D_%d",fManager->GetOutputEventNr());
281     fdigits->Write(nameDigits);
282     cout<<nameDigits<<endl;
283     wd->cd();
284   }
285 }
286
287
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 }