Digits and RecPoint with 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 #include "AliSTARTLoader.h"
41
42 #include <stdlib.h>
43 #include <Riostream.h>
44 #include <Riostream.h>
45
46 ClassImp(AliSTARTDigitizer)
47
48 //___________________________________________
49   AliSTARTDigitizer::AliSTARTDigitizer()  :AliDigitizer()
50 {
51 // Default ctor - don't use it
52   ;
53 }
54
55 //___________________________________________
56 AliSTARTDigitizer::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 //------------------------------------------------------------------------
69 AliSTARTDigitizer::~AliSTARTDigitizer()
70 {
71 // Destructor
72
73
74 }
75
76  //------------------------------------------------------------------------
77 Bool_t AliSTARTDigitizer::Init()
78 {
79 // Initialization
80  cout<<"AliSTARTDigitizer::Init"<<endl;
81  return kTRUE;
82 }
83  
84
85 //---------------------------------------------------------------------
86
87 void AliSTARTDigitizer::Exec(Option_t* option)
88 {
89
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");
96   cout<<" outgime "<<outgime<<endl;
97
98 #ifdef DEBUG
99   cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
100 #endif
101   //
102   // From hits to digits
103   //
104   Int_t hit, nhits;
105   Int_t CountEr[13],CountEl[13];                                                        //!!!
106   Int_t volume,pmt,tr,tl,sumRight;
107   Float_t timediff,timeav;
108   Float_t besttimeright,besttimeleft,meanTime;
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
115   //  Int_t thresholdAmpl=10;
116
117   ftimeRightTDC = new TArrayI(12); 
118   ftimeLeftTDC = new TArrayI(12); 
119   fRightADC = new TArrayI(12); 
120   fLeftADC = new TArrayI(12); 
121   
122   inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
123   inRL->LoadgAlice();
124   
125   //  fHits = new TClonesArray ("AliSTARThit", 1000);
126   fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000);                     //!!!
127   AliSTART *START  = (AliSTART*) gAlice->GetDetector("START");
128   AliSTARThit  *startHit;
129   AliSTARThitPhoton  *startHitPhoton;                                           //!!!
130   TBranch *brHits=0;
131   TBranch *brHitPhoton=0;
132   fdigits= new AliSTARTdigit();
133
134   Int_t nFiles=fManager->GetNinputs();
135   for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {
136
137     besttimeright=9999.;
138     besttimeleft=9999.;
139     Int_t timeDiff=0;
140     Int_t timeAv=0;
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       }
147
148     inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
149     //   AliSTARTLoader *pSTARTloader = (AliSTARTLoader*)fLoader;
150     //  pSTARTLoader->LoadHits("READ");
151     ingime = inRL->GetLoader("STARTLoader");
152     cout<<" ingime "<<ingime<<endl;
153     //    AliSTARTLoader *pSTARTLoader ;
154     ingime->LoadHits("READ");//probably it is necessary to load them before
155     ingime->LoadDigits("UPDATE");//probably it is necessary to load them before
156     TClonesArray *STARThitsPhotons = START->Photons ();
157     TClonesArray *fHits = START->Hits ();
158     //    cout<<" Load  "<<AliSTARTLoader::LoadDigits()<<endl;
159
160     TTree *th = ingime->TreeH();
161     brHits = th->GetBranch("START");
162     brHitPhoton = th->GetBranch("STARThitPhoton");
163     if (brHits) {
164       cout<<" brHits "<<endl;
165       START->SetHitsAddressBranch(brHits,brHitPhoton);
166     }else{
167       cerr<<"EXEC Branch START hit not found"<<endl;
168       exit(111);
169     } 
170     Int_t ntracks    = (Int_t) th->GetEntries();
171     cout<<" ntracks "<<ntracks<<endl;
172     if (ntracks<=0) return;
173     // Start loop on tracks in the photon hits containers 
174     // for amplitude
175     /*
176     if(brHitPhoton) {
177       cout<<"brHitPhoton "<<endl; 
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
194     */
195     // Start loop on tracks in the hits containers
196     cout<<" fHits "<<fHits<<endl;;
197     for (Int_t track=0; track<ntracks;track++) {
198       brHits->GetEntry(track);
199       nhits = fHits->GetEntriesFast();
200       //  cout<<" brHits hits "<<nhits<<endl;
201       for (hit=0;hit<nhits;hit++) {
202         startHit   = (AliSTARThit*) fHits->UncheckedAt(hit);
203         pmt=startHit->fPmt;
204         volume = startHit->fVolume;
205         if(volume==1){
206           timeright[pmt] = startHit->fTime;
207           if(timeright[pmt]<besttimeright)
208             //&&CountEr[pmt-1]>thresholdAmpl)
209             {
210             besttimeright=timeright[pmt];
211           } //timeright
212         }//time for right shoulder
213         if(volume==2){            
214           timeleft[pmt] = startHit->fTime;
215           if(timeleft[pmt]<besttimeleft)
216             //&&CountEl[pmt-1]>thresholdAmpl) 
217             {
218             besttimeleft=timeleft[pmt];
219             
220           } //timeleftbest
221         }//time for left shoulder
222       } //hit loop
223     } //track loop
224  
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
232     //folding with experimental time distribution
233     
234     besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
235     //    cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
236     bestRightADC=Int_t (besttimerightGaus*1000/channelWidth);
237     Float_t koef=69.7/350.;
238     besttimeleft=koef*besttimeleft;
239     besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
240     
241     bestLeftADC=Int_t (besttimeleftGaus*1000/channelWidth);
242     timediff=besttimerightGaus-besttimeleftGaus;
243     cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
244     meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
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       }
291     else
292       {timeAv=999999; timeDiff=99999;}
293
294 // trick to find out output dir:
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");
320   }
321 }
322
323
324 //------------------------------------------------------------------------
325 Bool_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 }
336 //----------------------------------------------------------------------------