1da215c4c60450c9115b49b9c8394b23649d616f
[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 <TFile.h>
20 #include <TDirectory.h>
21 #include <TRandom.h>
22 #include <TArrayI.h>
23 #include <TError.h>
24 #include <TH1F.h>
25
26 #include "AliLog.h"
27 #include "AliSTARTDigitizer.h"
28 #include "AliSTART.h"
29 #include "AliSTARThit.h"
30 #include "AliSTARTdigit.h"
31 #include "AliRunDigitizer.h"
32 #include <AliDetector.h>
33 #include "AliRun.h"
34 #include <AliLoader.h>
35 #include <AliRunLoader.h>
36 #include <stdlib.h>
37 #include <Riostream.h>
38 #include <Riostream.h>
39
40 ClassImp(AliSTARTDigitizer)
41
42 //___________________________________________
43   AliSTARTDigitizer::AliSTARTDigitizer()  :AliDigitizer()
44 {
45 // Default ctor - don't use it
46   ;
47 }
48
49 //___________________________________________
50 AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager) 
51   :AliDigitizer(manager),
52    fSTART(0),
53    fHits(0),
54    fdigits(0),
55    ftimeTDC(0),
56    fADC(0),
57    ftimeTDCAmp(0),
58    fADCAmp(0),
59    fSumMult(0),
60    fEff(0)
61 {
62   //    cout<<"AliSTARTDigitizer::AliSTARTDigitizer"<<endl;
63 // ctor which should be used
64
65   AliDebug(1,"processed");
66
67   fSTART = 0;
68   fHits = 0;
69   fdigits = 0;
70
71   ftimeTDC = new TArrayI(24); 
72   fADC = new TArrayI(24); 
73   ftimeTDCAmp = new TArrayI(24); 
74   fADCAmp = new TArrayI(24); 
75   fSumMult = new TArrayI(6); 
76   
77
78   TFile* file = TFile::Open("$ALICE_ROOT/START/PMTefficiency.root");
79   fEff = (TH1F*) file->Get("hEff")->Clone();
80   fEff->SetDirectory(NULL);
81   file->Close();
82   delete file;
83 }
84
85 //------------------------------------------------------------------------
86 AliSTARTDigitizer::~AliSTARTDigitizer()
87 {
88 // Destructor
89
90   AliDebug(1,"START");
91
92   delete ftimeTDC;
93   delete fADC;
94   delete fEff;
95   delete ftimeTDCAmp;
96   delete  fADCAmp;
97   delete fSumMult;
98 }
99
100  //------------------------------------------------------------------------
101 Bool_t AliSTARTDigitizer::Init()
102 {
103 // Initialization
104   AliDebug(1," Init");
105  return kTRUE;
106 }
107  
108
109 //---------------------------------------------------------------------
110
111 void AliSTARTDigitizer::Exec(Option_t* /*option*/)
112 {
113
114   /*
115     Produde digits from hits
116         digits is TObject and includes
117         We are writing array if left & right  TDC
118         left & right  ADC (will need for slow simulation)
119         TOF first particle left & right
120         mean time and time difference (vertex position)
121         
122   */
123
124   //output loader 
125   AliRunLoader *outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
126   AliLoader * pOutStartLoader = outRL->GetLoader("STARTLoader");
127
128   AliDebug(1,"start...");
129   cout<<" AliSTARTDigitizer::Exec "<<endl;
130   //input loader
131   //
132   // From hits to digits
133   //
134   Int_t hit, nhits;
135   Int_t countE[24], sumMult[3];
136   Int_t volume,pmt,tr,qt,qtAmp;
137   Int_t  bestRightTDC,bestLeftTDC;
138   Float_t time[24],besttime[24], timeGaus[24] ;
139   Float_t channelWidth=25.; //ps
140     //Q->T-> coefficients !!!! should be asked!!!
141   Float_t ph2mV = 150./500.;
142   Int_t threshold =50; //photoelectrons
143   Int_t mV2channel=200000/(25*25);  //5V -> 200ns
144
145   
146   AliSTARThit  *startHit;
147   TBranch *brHits=0;
148
149   Int_t nFiles=fManager->GetNinputs();
150   for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {
151     if (inputFile < nFiles-1) {
152       AliWarning(Form("ignoring input stream %d", inputFile));
153       continue;
154       
155     }
156     
157     Float_t besttimeright=99999.;
158     Float_t besttimeleft=99999.;
159     Int_t timeDiff, meanTime;
160     
161     for (Int_t i0=0; i0<24; i0++)
162       {
163         time[i0]=99999;  besttime[i0]=99999;    countE[i0]=0;
164       }
165     for ( Int_t imu=0; imu<3; imu++) sumMult[imu]=0;
166     AliRunLoader * inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
167     AliLoader * pInStartLoader = inRL->GetLoader("STARTLoader");
168     if (!inRL->GetAliRun()) inRL->LoadgAlice();
169     AliSTART *fSTART  = (AliSTART*)inRL ->GetAliRun()->GetDetector("START");
170
171        //read Hits 
172     pInStartLoader->LoadHits("READ");//probably it is necessary to load them before
173     TClonesArray *fHits = fSTART->Hits ();
174     TTree *th = pInStartLoader->TreeH();
175     brHits = th->GetBranch("START");
176     if (brHits) {
177       fSTART->SetHitsAddressBranch(brHits);
178     }else{
179       AliError("Branch START hit not found");
180       exit(111);
181     } 
182     Int_t ntracks    = (Int_t) th->GetEntries();
183     
184     if (ntracks<=0) return;
185     // Start loop on tracks in the hits containers
186     for (Int_t track=0; track<ntracks;track++) {
187       brHits->GetEntry(track);
188       nhits = fHits->GetEntriesFast();
189       for (hit=0;hit<nhits;hit++) 
190         {
191           startHit   = (AliSTARThit*) fHits->UncheckedAt(hit);
192           if (!startHit) {
193             AliError("The unchecked hit doesn't exist");
194             break;
195           }
196           pmt=startHit->Pmt();
197           Int_t numpmt=pmt-1;
198                   Double_t e=startHit->Etot();
199           volume = startHit->Volume();
200           
201           if(e>0 && RegisterPhotoE(e)) {
202             countE[numpmt]++;
203             besttime[numpmt] = startHit->Time();
204             if(besttime[numpmt]<time[numpmt])
205               {
206                 time[numpmt]=besttime[numpmt];
207               }
208           } 
209         } //hits loop
210     } //track loop
211     
212     //spread time right&left by 25ps   && besttime
213     for (Int_t ipmt=0; ipmt<12; ipmt++){
214       if(countE[ipmt] > threshold) {
215         timeGaus[ipmt]=gRandom->Gaus(time[ipmt],25);
216         if(timeGaus[ipmt]<besttimeleft) besttimeleft=timeGaus[ipmt]; //timeleft
217         sumMult[0] +=  countE[ipmt];
218         sumMult[1] +=  countE[ipmt];
219       }
220     }
221     for ( Int_t ipmt=12; ipmt<24; ipmt++){
222       if(countE[ipmt] > threshold) {
223         timeGaus[ipmt]=gRandom->Gaus(time[ipmt],25); 
224         if(timeGaus[ipmt]<besttimeright)  besttimeright=timeGaus[ipmt]; //timeright
225         sumMult[0] +=  countE[ipmt];
226         sumMult[2] +=  countE[ipmt];
227       } 
228     }
229     //ADC features fill digits
230    //folding with experimental time distribution
231     Float_t c = 0.0299792; // cm/ps
232     Float_t koef=(350.-69.7)/c;
233     Float_t  besttimeleftR= besttimeleft;
234     besttimeleft=koef+besttimeleft;
235     bestLeftTDC=Int_t (besttimeleftR/channelWidth);
236     bestRightTDC=Int_t (besttimeright/channelWidth);
237     timeDiff=Int_t ((besttimeright-besttimeleftR)/channelWidth);
238     meanTime=Int_t (((besttimeright+besttimeleft)/2.)/channelWidth);
239     AliDebug(2,Form(" time in ns %f ", besttimeleft));
240     for (Int_t i=0; i<24; i++)
241       {
242         //  fill TDC
243         tr= Int_t (timeGaus[i]/channelWidth); 
244         if(timeGaus[i]>50000) tr=0;
245         
246         //fill ADC
247         Int_t al= countE[i]; 
248         // QTC procedure:
249         // phe -> mV 0.3; 1MIP ->500phe -> ln (amp (mV)) = 5;
250         // max 200ns, HIJING  mean 50000phe -> 15000mv -> ln = 15 (s zapasom)
251         // channel 25ps
252         if (al>threshold) {
253           qt=Int_t (TMath::Log(al*ph2mV) * mV2channel); 
254           qtAmp=Int_t (TMath::Log(al*10*ph2mV) * mV2channel);
255           cout<<i<<" "<<qt<<" "<<qtAmp<<endl;
256           fADC->AddAt(qt,i);
257           ftimeTDC->AddAt(tr,i);
258           fADCAmp->AddAt(qtAmp,i);
259           ftimeTDCAmp->AddAt(tr,i); //now is the same as non-amplified but will be change
260         }
261       } //pmt loop
262     for (Int_t im=0; im<3; im++)
263       { 
264         if (sumMult[im]>threshold){
265           Int_t sum=Int_t (TMath::Log(sumMult[im]*ph2mV) * mV2channel);
266           Int_t sumAmp=Int_t (TMath::Log(10*sumMult[im]*ph2mV) * mV2channel);
267           fSumMult->AddAt(sum,im);
268           fSumMult->AddAt(sumAmp,im+3);
269         }
270       } 
271
272     fSTART->AddDigit(bestRightTDC,bestLeftTDC,meanTime,timeDiff,fSumMult,
273                      ftimeTDC,fADC,ftimeTDCAmp,fADCAmp);
274     pOutStartLoader->UnloadHits();
275   } //input streams loop
276   
277     //load digits    
278     pOutStartLoader->LoadDigits("UPDATE");
279     TTree *treeD  = pOutStartLoader->TreeD();
280     if (treeD == 0x0) {
281       pOutStartLoader->MakeTree("D");
282       treeD = pOutStartLoader->TreeD();
283     }
284     AliSTART *fSTART  = (AliSTART*)outRL ->GetAliRun()->GetDetector("START");
285     fSTART->MakeBranch("D");
286      treeD->Reset();
287      treeD->Fill();
288   
289   pOutStartLoader->WriteDigits("OVERWRITE");
290   
291   fSTART->ResetDigits();
292   pOutStartLoader->UnloadDigits();
293     
294 }
295
296
297 //------------------------------------------------------------------------
298 Bool_t AliSTARTDigitizer::RegisterPhotoE(Double_t energy)
299 {
300
301   
302   //  Float_t hc=197.326960*1.e6; //mev*nm
303   Double_t hc=1.973*1.e-6; //gev*nm
304   //  cout<<"AliSTARTDigitizer::RegisterPhotoE >> energy "<<energy<<endl;
305   Float_t lambda=hc/energy;
306   Int_t bin=  fEff->GetXaxis()->FindBin(lambda);
307   Float_t eff=fEff->GetBinContent(bin);
308   Double_t  p = gRandom->Rndm();
309   if (p > eff)
310     return kFALSE;
311   
312   return kTRUE;
313 }
314 //----------------------------------------------------------------------------