Transition to NewIO
[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 #include "AliLoader.h"
39 #include "AliRunLoader.h"
40
41 #include <stdlib.h>
42 #include <Riostream.h>
43 #include <Riostream.h>
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
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
96 #ifdef DEBUG
97   cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
98 #endif
99   //
100   // From hits to digits
101   //
102   Int_t hit, nhits;
103   Int_t CountEr[13],CountEl[13];                                                        //!!!
104   Int_t volume,pmt,tr,tl,sumRight;
105   Float_t timediff,timeav;
106   Float_t besttimeright,besttimeleft,meanTime;
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
113   //  Int_t thresholdAmpl=10;
114
115   ftimeRightTDC = new TArrayI(12); 
116   ftimeLeftTDC = new TArrayI(12); 
117   fRightADC = new TArrayI(12); 
118   fLeftADC = new TArrayI(12); 
119   
120
121   fHits = new TClonesArray ("AliSTARThit", 1000);
122   fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000);                     //!!!
123   AliSTART *START  = (AliSTART*) gAlice->GetDetector("START");
124
125   AliSTARThit  *startHit;
126   AliSTARThitPhoton  *startHitPhoton;                                           //!!!
127   TBranch *brHits=0;
128   TBranch *brHitPhoton=0;
129   fdigits= new AliSTARTdigit();
130
131   Int_t nFiles=fManager->GetNinputs();
132   for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {
133
134  
135
136     besttimeright=9999.;
137     besttimeleft=9999.;
138     Int_t timeDiff=0;
139     Int_t timeAv=0;
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       }
146     TClonesArray *STARThits = START->Hits ();
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
151     TClonesArray *STARThitsPhotons = START->Photons ();
152
153
154    TTree *th = ingime->TreeH();
155     brHits = th->GetBranch("START");
156     brHitPhoton = th->GetBranch("STARThitPhoton");
157     if (brHits) {
158       START->SetHitsAddressBranch(brHits,brHitPhoton);
159     }else{
160       cerr<<"EXEC Branch START hit not found"<<endl;
161       exit(111);
162     } 
163     Int_t ntracks    = (Int_t) th->GetEntries();
164     cout<<" ntracks "<<ntracks<<endl;
165     if (ntracks<=0) return;
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
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){
195           timeright[pmt] = startHit->fTime;
196           if(timeright[pmt]<besttimeright)
197             //&&CountEr[pmt-1]>thresholdAmpl)
198             {
199             besttimeright=timeright[pmt];
200           } //timeright
201         }//time for right shoulder
202         if(volume==2){            
203           timeleft[pmt] = startHit->fTime;
204           if(timeleft[pmt]<besttimeleft)
205             //&&CountEl[pmt-1]>thresholdAmpl) 
206             {
207             besttimeleft=timeleft[pmt];
208             
209           } //timeleftbest
210         }//time for left shoulder
211       } //hit loop
212     } //track loop
213  
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
221     //folding with experimental time distribution
222     
223     besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
224     //    cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
225     bestRightADC=Int_t (besttimerightGaus*1000/channelWidth);
226     Float_t koef=69.7/350.;
227     besttimeleft=koef*besttimeleft;
228     besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
229     
230     bestLeftADC=Int_t (besttimeleftGaus*1000/channelWidth);
231     timediff=besttimerightGaus-besttimeleftGaus;
232     cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
233     meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
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       }
280     else
281       {timeAv=999999; timeDiff=99999;}
282
283 // trick to find out output dir:
284     outgime->WriteDigits("OVERWRITE");
285   }
286 }
287
288
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 }