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