New AliTRDcluster constructor
[u/mrichter/AliRoot.git] / TOF / AliTOFSDigitizer.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //__________________________________________________________//
19 //                                                          //
20 //   This is a TTask that constructs SDigits out of Hits    //
21 //   A Summable Digits is the "sum" of all hits in a pad    //
22 //   Detector response has been simulated via the method    //
23 //   SimulateDetectorResponse                               //
24 //                                                          //
25 //  -- Authors: F. Pierella, A. De Caro                     //
26 //   Use case: see AliTOFhits2sdigits.C macro in the CVS    //
27 //__________________________________________________________//
28
29 #include <TBenchmark.h>
30 #include <TClonesArray.h>
31 #include <TF1.h>
32 #include <TFile.h>
33 #include <TParticle.h>
34 #include <TTree.h>
35 #include <TRandom.h>
36 #include <TROOT.h>
37
38 #include "AliLoader.h"
39 #include "AliLog.h"
40 #include "AliMC.h"
41 #include "AliRunLoader.h"
42 #include "AliRun.h"
43
44 #include "AliTOFGeometry.h"
45 #include "AliTOFHitMap.h"
46 #include "AliTOFhitT0.h"
47 #include "AliTOFhit.h"
48 #include "AliTOFSDigitizer.h"
49 #include "AliTOFSDigit.h"
50 #include "AliTOF.h"
51
52
53 ClassImp(AliTOFSDigitizer)
54
55 //____________________________________________________________________________ 
56 AliTOFSDigitizer::AliTOFSDigitizer():
57   TTask("TOFSDigitizer",""),
58   fEvent1(-1),
59   fEvent2(-1),
60   ftail(0x0),
61   fHeadersFile(""),
62   fRunLoader(0x0),
63   fTOFLoader(0x0),
64   fTOFGeometry(new AliTOFGeometry()),
65   fSelectedSector(-1), 
66   fSelectedPlate(-1),
67   fTimeResolution(0),
68   fpadefficiency(0),
69   fEdgeEffect(-1),
70   fEdgeTails(-1),
71   fHparameter(0),
72   fH2parameter(0),
73   fKparameter(0),
74   fK2parameter(0),
75   fEffCenter(0),
76   fEffBoundary(0),
77   fEff2Boundary(0),
78   fEff3Boundary(0),
79   fAddTRes(0),
80   fResCenter(0),
81   fResBoundary(0),
82   fResSlope(0),
83   fTimeWalkCenter(0),
84   fTimeWalkBoundary(0),
85   fTimeWalkSlope(0),
86   fTimeDelayFlag(-1),
87   fPulseHeightSlope(0),
88   fTimeDelaySlope(0),
89   fMinimumCharge(0),
90   fChargeSmearing(0),
91   fLogChargeSmearing(0),
92   fTimeSmearing(0),
93   fAverageTimeFlag(-1),
94   fAdcBin(0),
95   fAdcMean(0),
96   fAdcRms(0)
97 {
98   // ctor
99 }
100
101 //------------------------------------------------------------------------
102 AliTOFSDigitizer::AliTOFSDigitizer(const AliTOFSDigitizer &source):
103   TTask(source),
104   fEvent1(-1),
105   fEvent2(-1),
106   ftail(0x0),
107   fHeadersFile(""),
108   fRunLoader(0x0),
109   fTOFLoader(0x0),
110   fTOFGeometry(0x0),
111   fSelectedSector(-1), 
112   fSelectedPlate(-1),
113   fTimeResolution(0),
114   fpadefficiency(0),
115   fEdgeEffect(-1),
116   fEdgeTails(-1),
117   fHparameter(0),
118   fH2parameter(0),
119   fKparameter(0),
120   fK2parameter(0),
121   fEffCenter(0),
122   fEffBoundary(0),
123   fEff2Boundary(0),
124   fEff3Boundary(0),
125   fAddTRes(0),
126   fResCenter(0),
127   fResBoundary(0),
128   fResSlope(0),
129   fTimeWalkCenter(0),
130   fTimeWalkBoundary(0),
131   fTimeWalkSlope(0),
132   fTimeDelayFlag(-1),
133   fPulseHeightSlope(0),
134   fTimeDelaySlope(0),
135   fMinimumCharge(0),
136   fChargeSmearing(0),
137   fLogChargeSmearing(0),
138   fTimeSmearing(0),
139   fAverageTimeFlag(-1),
140   fAdcBin(0),
141   fAdcMean(0),
142   fAdcRms(0)
143 {
144   // copy constructor
145   this->fTOFGeometry=source.fTOFGeometry;
146
147 }
148
149 //____________________________________________________________________________ 
150 AliTOFSDigitizer& AliTOFSDigitizer::operator=(const AliTOFSDigitizer &source)
151 {
152   // ass. op.
153   this->fTOFGeometry=source.fTOFGeometry;
154   return *this;
155
156 }
157
158 //____________________________________________________________________________ 
159 AliTOFSDigitizer::AliTOFSDigitizer(const char* HeaderFile, Int_t evNumber1, Int_t nEvents):
160   TTask("TOFSDigitizer",""),
161   fEvent1(-1),
162   fEvent2(-1),
163   ftail(0x0),
164   fHeadersFile(HeaderFile), // input filename (with hits)
165   fRunLoader(0x0),
166   fTOFLoader(0x0),
167   fTOFGeometry(0x0),
168   fSelectedSector(-1), // by default we sdigitize all sectors
169   fSelectedPlate(-1),  // by default we sdigitize all plates in all sectors
170   fTimeResolution(0),
171   fpadefficiency(0),
172   fEdgeEffect(-1),
173   fEdgeTails(-1),
174   fHparameter(0),
175   fH2parameter(0),
176   fKparameter(0),
177   fK2parameter(0),
178   fEffCenter(0),
179   fEffBoundary(0),
180   fEff2Boundary(0),
181   fEff3Boundary(0),
182   fAddTRes(0),
183   fResCenter(0),
184   fResBoundary(0),
185   fResSlope(0),
186   fTimeWalkCenter(0),
187   fTimeWalkBoundary(0),
188   fTimeWalkSlope(0),
189   fTimeDelayFlag(-1),
190   fPulseHeightSlope(0),
191   fTimeDelaySlope(0),
192   fMinimumCharge(0),
193   fChargeSmearing(0),
194   fLogChargeSmearing(0),
195   fTimeSmearing(0),
196   fAverageTimeFlag(-1),
197   fAdcBin(0),
198   fAdcMean(0),
199   fAdcRms(0)
200 {
201   //ctor, reading from input file 
202   
203   TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data());
204   
205   //File was not opened yet open file and get alirun object
206   if (file == 0) {
207     file   = TFile::Open(fHeadersFile.Data(),"update") ;
208     gAlice = (AliRun *) file->Get("gAlice") ;
209   }
210   
211   // add Task to //root/Tasks folder
212   TString evfoldname = AliConfig::GetDefaultEventFolderName();
213   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
214   if (!fRunLoader)
215     fRunLoader = AliRunLoader::Open(HeaderFile);//open session and mount on default event folder
216   if (fRunLoader == 0x0)
217     {
218       AliFatal("Event is not loaded. Exiting");
219       return;
220     }
221
222   fRunLoader->CdGAFile();
223   TDirectory *savedir=gDirectory;
224   TFile *in=(TFile*)gFile;
225
226   if (!in->IsOpen()) {
227     AliWarning("Geometry file is not open default TOF geometry will be used");
228     fTOFGeometry = new AliTOFGeometry();
229   }
230   else {
231     in->cd();
232     fTOFGeometry = (AliTOFGeometry*)in->Get("TOFgeometry");
233   }
234
235   savedir->cd();
236
237   if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader();
238   
239   if (evNumber1>=0) fEvent1 = evNumber1;
240   else fEvent1=0;
241   
242   if (nEvents==0) fEvent2 = (Int_t)(fRunLoader->GetNumberOfEvents());
243   else if (nEvents>0) fEvent2 = evNumber1+nEvents;
244   else fEvent2 = 1;
245   
246   if (!(fEvent2>fEvent1)) {
247     AliError(Form("fEvent2 = %d <= fEvent1 = %d", fEvent2, fEvent1));
248     fEvent1 = 0;
249     fEvent2 = 1;
250     AliError(Form("Correction: fEvent2 = %d <= fEvent1 = %d", fEvent2, fEvent1));
251   }
252   
253   // init parameters for sdigitization
254   InitParameters();
255   
256   fTOFLoader = fRunLoader->GetLoader("TOFLoader");
257   if (fTOFLoader == 0x0)
258     {
259       AliFatal("Can not find TOF loader in event. Exiting.");
260       return;
261     }
262   fTOFLoader->PostSDigitizer(this);
263 }
264
265 //____________________________________________________________________________ 
266 AliTOFSDigitizer::~AliTOFSDigitizer()
267 {
268   // dtor
269   fTOFLoader->CleanSDigitizer();
270
271   delete fTOFGeometry;
272
273 }
274
275 //____________________________________________________________________________ 
276 void AliTOFSDigitizer::InitParameters()
277 {
278   // set parameters for detector simulation
279   
280   fTimeResolution = 0.080; //0.120; OLD
281   fpadefficiency  = 0.99 ;
282   fEdgeEffect     = 2   ;
283   fEdgeTails      = 0   ;
284   fHparameter     = 0.4 ;
285   fH2parameter    = 0.15;
286   fKparameter     = 0.5 ;
287   fK2parameter    = 0.35;
288   fEffCenter      = fpadefficiency;
289   fEffBoundary    = 0.65;
290   fEff2Boundary   = 0.90;
291   fEff3Boundary   = 0.08;
292   fAddTRes        = 68. ; // \sqrt{2x20^2 + 15^2 + 2x10^2 + 30^2 + 50^2} (p-p)
293   //fAddTRes      = 48. ; // \sqrt{2x20^2 + 15^2 + 2x10^2 + 30^2 + 15^2} (Pb-Pb)
294   // 30^2+20^2+40^2+50^2+50^2+50^2 = 10400 ps^2 (very old value)
295   fResCenter      = 35. ; //50. ; // OLD
296   fResBoundary    = 70. ;
297   fResSlope       = 37. ; //40. ; // OLD
298   fTimeWalkCenter = 0.  ;
299   fTimeWalkBoundary=0.  ;
300   fTimeWalkSlope  = 0.  ;
301   fTimeDelayFlag  = 1   ;
302   fPulseHeightSlope=2.0 ;
303   fTimeDelaySlope =0.060;
304   // was fMinimumCharge = TMath::Exp(fPulseHeightSlope*fKparameter/2.);
305   fMinimumCharge = TMath::Exp(-fPulseHeightSlope*fHparameter);
306   fChargeSmearing=0.0   ;
307   fLogChargeSmearing=0.13;
308   fTimeSmearing   =0.022;
309   fAverageTimeFlag=0    ;
310
311   fAdcBin   = 0.25;    // 1 ADC bin = 0.25 pC (or 0.03 pC)
312   fAdcMean  = 50.;     // ADC distribution mpv value for Landau (in bins)
313                        // it corresponds to a mean value of ~100 bins
314   fAdcRms   = 25.;     // ADC distribution rms value (in bins)
315                        // it corresponds to distribution rms ~50 bins
316 }
317
318 //__________________________________________________________________
319 Double_t TimeWithTail(Double_t* x, Double_t* par)
320 {
321   // sigma - par[0], alpha - par[1], part - par[2]
322   //  at x<part*sigma - gauss
323   //  at x>part*sigma - TMath::Exp(-x/alpha)
324   Float_t xx =x[0];
325   Double_t f;
326   if(xx<par[0]*par[2]) {
327     f = TMath::Exp(-xx*xx/(2*par[0]*par[0]));
328   } else {
329     f = TMath::Exp(-(xx-par[0]*par[2])/par[1]-0.5*par[2]*par[2]);
330   }
331   return f;
332 }
333
334 //____________________________________________________________________________
335 void AliTOFSDigitizer::Exec(Option_t *verboseOption) { 
336   //execute TOF sdigitization
337   if (strstr(verboseOption,"tim") || strstr(verboseOption,"all"))
338     gBenchmark->Start("TOFSDigitizer");
339
340   if (fEdgeTails) ftail = new TF1("tail",TimeWithTail,-2,2,3);
341   
342   Int_t nselectedHits=0;
343   Int_t ntotalsdigits=0;
344   Int_t ntotalupdates=0;
345   Int_t nnoisesdigits=0;
346   Int_t nsignalsdigits=0;
347   Int_t nHitsFromPrim=0;
348   Int_t nHitsFromSec=0;
349   Int_t nlargeTofDiff=0;
350
351   Bool_t thereIsNotASelection=(fSelectedSector==-1) && (fSelectedPlate==-1);
352
353   if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
354   gAlice = fRunLoader->GetAliRun();
355
356   fRunLoader->LoadKinematics();
357   
358   AliTOF *tof = (AliTOF *) gAlice->GetDetector("TOF");
359   
360   if (!tof) {
361     AliError("TOF not found");
362     return;
363   }
364   
365   fTOFLoader->LoadHits("read");
366   fTOFLoader->LoadSDigits("recreate");
367   
368   for (Int_t iEvent=fEvent1; iEvent<fEvent2; iEvent++) {
369     //AliInfo(Form("------------------- %s -------------", GetName()));
370     //AliInfo(Form("Sdigitizing event %i", iEvent));
371
372     fRunLoader->GetEvent(iEvent);
373
374     TTree *hitTree = fTOFLoader->TreeH ();
375     if (!hitTree) return;
376
377     if (fTOFLoader->TreeS () == 0) fTOFLoader->MakeTree ("S");
378     
379     //Make branch for digits
380     tof->MakeBranch("S");
381     
382     // recreate TClonesArray fSDigits - for backward compatibility
383     if (tof->SDigits() == 0) {
384       tof->CreateSDigitsArray();
385     } else {
386       tof->RecreateSDigitsArray();
387     }
388
389     tof->SetTreeAddress();
390
391     Int_t version=tof->IsVersion();
392
393     Int_t nselectedHitsinEv=0;
394     Int_t ntotalsdigitsinEv=0;
395     Int_t ntotalupdatesinEv=0;
396     Int_t nnoisesdigitsinEv=0;
397     Int_t nsignalsdigitsinEv=0;
398
399     TParticle *particle;
400     //AliTOFhit *tofHit;
401     TClonesArray *tofHitArray = tof->Hits();
402
403     // create hit map
404     AliTOFHitMap *hitMap = new AliTOFHitMap(tof->SDigits(), fTOFGeometry);
405
406     TBranch * tofHitsBranch = hitTree->GetBranch("TOF");
407
408     Int_t ntracks = static_cast<Int_t>(hitTree->GetEntries());
409     for (Int_t track = 0; track < ntracks; track++)
410     {
411       gAlice->ResetHits();
412       tofHitsBranch->GetEvent(track);
413
414       AliMC *mcApplication = (AliMC*)gAlice->GetMCApp();
415
416       particle = mcApplication->Particle(track);
417       Int_t nhits = tofHitArray->GetEntriesFast();
418       // cleaning all hits of the same track in the same pad volume
419       // it is a rare event, however it happens
420
421       Int_t previousTrack =-1;
422       Int_t previousSector=-1;
423       Int_t previousPlate =-1;
424       Int_t previousStrip =-1;
425       Int_t previousPadX  =-1;
426       Int_t previousPadZ  =-1;
427
428       for (Int_t hit = 0; hit < nhits; hit++) {
429         Int_t    vol[5];       // location for a digit
430         Int_t  digit[2];     // TOF digit variables
431         Int_t tracknum;
432         Float_t dxPad;
433         Float_t dzPad;
434         Float_t geantTime;
435
436         // fp: really sorry for this, it is a temporary trick to have
437         // track length too
438         if (version<6) { //(version!=6 && version!=7)
439           AliTOFhit *tofHit = (AliTOFhit *) tofHitArray->UncheckedAt(hit);
440           tracknum = tofHit->GetTrack();
441           vol[0] = tofHit->GetSector();
442           vol[1] = tofHit->GetPlate();
443           vol[2] = tofHit->GetStrip();
444           vol[3] = tofHit->GetPadx();
445           vol[4] = tofHit->GetPadz();
446           dxPad = tofHit->GetDx();
447           dzPad = tofHit->GetDz();
448           geantTime = tofHit->GetTof(); // unit [s]
449         } else {
450           AliTOFhitT0 *tofHit = (AliTOFhitT0 *) tofHitArray->UncheckedAt(hit);
451           tracknum = tofHit->GetTrack();
452           vol[0] = tofHit->GetSector();
453           vol[1] = tofHit->GetPlate();
454           vol[2] = tofHit->GetStrip();
455           vol[3] = tofHit->GetPadx();
456           vol[4] = tofHit->GetPadz();
457           dxPad = tofHit->GetDx();
458           dzPad = tofHit->GetDz();
459           geantTime = tofHit->GetTof(); // unit [s]
460         }
461         
462         geantTime *= 1.e+09;  // conversion from [s] to [ns]
463         
464         // selection case for sdigitizing only hits in a given plate of a given sector
465         if(thereIsNotASelection || (vol[0]==fSelectedSector && vol[1]==fSelectedPlate)){
466           
467           Bool_t dummy=((tracknum==previousTrack) && (vol[0]==previousSector) && (vol[1]==previousPlate) && (vol[2]==previousStrip));
468           
469           Bool_t isCloneOfThePrevious=dummy && ((vol[3]==previousPadX) && (vol[4]==previousPadZ));
470           
471           Bool_t isNeighOfThePrevious=dummy && ((((vol[3]==previousPadX-1) || (vol[3]==previousPadX+1)) && (vol[4]==previousPadZ)) || ((vol[3]==previousPadX) && ((vol[4]==previousPadZ+1) || (vol[4]==previousPadZ-1))));
472           
473           if(!isCloneOfThePrevious && !isNeighOfThePrevious){
474             // update "previous" values
475             // in fact, we are yet in the future, so the present is past
476             previousTrack=tracknum;
477             previousSector=vol[0];
478             previousPlate=vol[1];
479             previousStrip=vol[2];
480             previousPadX=vol[3];
481             previousPadZ=vol[4];
482             
483             nselectedHits++;
484             nselectedHitsinEv++;
485             if (particle->GetFirstMother() < 0) nHitsFromPrim++; // counts hits due to primary particles
486             
487             Float_t xStrip=AliTOFGeometry::XPad()*(vol[3]+0.5-0.5*AliTOFGeometry::NpadX())+dxPad;
488             Float_t zStrip=AliTOFGeometry::ZPad()*(vol[4]+0.5-0.5*AliTOFGeometry::NpadZ())+dzPad;
489
490             Int_t nActivatedPads = 0, nFiredPads = 0;
491             Bool_t isFired[4] = {kFALSE, kFALSE, kFALSE, kFALSE};
492             Float_t tofAfterSimul[4] = {0., 0., 0., 0.};
493             Float_t qInduced[4] = {0.,0.,0.,0.};
494             Int_t nPlace[4] = {0, 0, 0, 0};
495             Float_t averageTime = 0.;
496             SimulateDetectorResponse(zStrip,xStrip,geantTime,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime);
497             if(nFiredPads) {
498               for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
499                 if(isFired[indexOfPad]){ // the pad has fired
500                   Float_t timediff=geantTime-tofAfterSimul[indexOfPad];
501                   
502                   if(timediff>=0.2) nlargeTofDiff++;
503                   
504                   digit[0] = (Int_t) ((tofAfterSimul[indexOfPad]*1.e+03)/AliTOFGeometry::TdcBinWidth()); // TDC bin number (each bin -> 24.4 ps)
505                   
506                   Float_t landauFactor = gRandom->Landau(fAdcMean, fAdcRms); 
507                   digit[1] = (Int_t) (qInduced[indexOfPad] * landauFactor); // ADC bins (each bin -> 0.25 (or 0.03) pC)
508
509                   // recalculate the volume only for neighbouring pads
510                   if(indexOfPad){
511                     (nPlace[indexOfPad]<=AliTOFGeometry::NpadX()) ? vol[4] = 0 : vol[4] = 1;
512                     (nPlace[indexOfPad]<=AliTOFGeometry::NpadX()) ? vol[3] = nPlace[indexOfPad] - 1 : vol[3] = nPlace[indexOfPad] - AliTOFGeometry::NpadX() - 1;
513                   }
514                   // check if two sdigit are on the same pad;
515                   // in that case we sum the two or more sdigits
516                   if (hitMap->TestHit(vol) != kEmpty) {
517                     AliTOFSDigit *sdig = static_cast<AliTOFSDigit*>(hitMap->GetHit(vol));
518                     Int_t tdctime = (Int_t) digit[0];
519                     Int_t adccharge = (Int_t) digit[1];
520                     sdig->Update(AliTOFGeometry::TdcBinWidth(),tdctime,adccharge,tracknum);
521                     ntotalupdatesinEv++;
522                     ntotalupdates++;
523                   } else {
524                     
525                     tof->AddSDigit(tracknum, vol, digit);
526                     
527                     if(indexOfPad){
528                       nnoisesdigits++;
529                       nnoisesdigitsinEv++;
530                     } else {
531                       nsignalsdigits++;
532                       nsignalsdigitsinEv++;
533                     }
534                     ntotalsdigitsinEv++;  
535                     ntotalsdigits++;
536                     hitMap->SetHit(vol);
537                   } // if (hitMap->TestHit(vol) != kEmpty)
538                 } // if(isFired[indexOfPad])
539               } // end loop on nActivatedPads
540             } // if(nFiredPads) i.e. if some pads has fired
541           } // close if(!isCloneOfThePrevious)
542         } // close the selection on sector and plate
543       } // end loop on hits for the current track
544     } // end loop on ntracks
545     
546     delete hitMap;
547     
548     fTOFLoader->TreeS()->Reset();
549     fTOFLoader->TreeS()->Fill();
550     fTOFLoader->WriteSDigits("OVERWRITE");
551     
552     if (tof->SDigits()) tof->ResetSDigits();
553     
554     if (strstr(verboseOption,"all")) {
555       AliInfo("----------------------------------------");
556       AliInfo("       <AliTOFSDigitizer>    ");
557       AliInfo(Form("After sdigitizing %d hits in event %d", nselectedHitsinEv, iEvent));
558       //" (" << nHitsFromPrim << " from primaries and " << nHitsFromSec << " from secondaries) TOF hits, " 
559       AliInfo(Form("%d digits have been created", ntotalsdigitsinEv));
560       AliInfo(Form("(%d due to signals and %d due to border effect)", nsignalsdigitsinEv, nnoisesdigitsinEv));
561       AliInfo(Form("%d total updates of the hit map have been performed in current event", ntotalupdatesinEv));
562       AliInfo("----------------------------------------");
563     }
564
565   } //event loop on events
566
567     fTOFLoader->UnloadSDigits();
568     fTOFLoader->UnloadHits();
569     fRunLoader->UnloadKinematics();
570     //fRunLoader->UnloadgAlice();
571
572   // free used memory
573   if (ftail){
574     delete ftail;
575     ftail = 0;
576   }
577   
578   nHitsFromSec=nselectedHits-nHitsFromPrim;
579   if(strstr(verboseOption,"all")){
580     AliInfo("----------------------------------------");
581     AliInfo("----------------------------------------");
582     AliInfo("-----------SDigitization Summary--------");
583     AliInfo("       <AliTOFSDigitizer>     ");
584     AliInfo(Form("After sdigitizing %d hits", nselectedHits));
585     AliInfo(Form("in %d events", fEvent2-fEvent1));
586 //" (" << nHitsFromPrim << " from primaries and " << nHitsFromSec << " from secondaries) TOF hits, " 
587     AliInfo(Form("%d sdigits have been created", ntotalsdigits));
588     AliInfo(Form("(%d due to signals and " 
589                  "%d due to border effect)", nsignalsdigits, nnoisesdigits));
590     AliInfo(Form("%d total updates of the hit map have been performed", ntotalupdates));
591     AliInfo(Form("in %d cases the time of flight difference is greater than 200 ps", nlargeTofDiff));
592   }
593
594
595   if(strstr(verboseOption,"tim") || strstr(verboseOption,"all")){
596     gBenchmark->Stop("TOFSDigitizer");
597     AliInfo("AliTOFSDigitizer:");
598     AliInfo(Form("   took %f seconds in order to make sdigits " 
599          "%f seconds per event", gBenchmark->GetCpuTime("TOFSDigitizer"), gBenchmark->GetCpuTime("TOFSDigitizer")/(fEvent2-fEvent1)));
600     AliInfo(" +++++++++++++++++++++++++++++++++++++++++++++++++++ ");
601   }
602
603 }
604
605 //__________________________________________________________________
606 void AliTOFSDigitizer::Print(Option_t* /*opt*/)const
607 {
608   AliInfo(Form(" ------------------- %s ------------- ", GetName()));
609 }
610
611 //__________________________________________________________________
612 void AliTOFSDigitizer::SelectSectorAndPlate(Int_t sector, Int_t plate)
613 {
614   //Select sector and plate
615   Bool_t isaWrongSelection=(sector < 0) || (sector >= AliTOFGeometry::NSectors()) || (plate < 0) || (plate >= AliTOFGeometry::NPlates());
616   if(isaWrongSelection){
617     AliError("You have selected an invalid value for sector or plate ");
618     AliError(Form("The correct range for sector is [0,%d]", AliTOFGeometry::NSectors()-1));
619     AliError(Form("The correct range for plate  is [0,%d]",  AliTOFGeometry::NPlates()-1));
620     AliError("By default we continue sdigitizing all hits in all plates of all sectors");
621   } else {
622     fSelectedSector=sector;
623     fSelectedPlate =plate;
624     AliInfo(Form("SDigitizing only hits in plate %d of the sector %d", fSelectedPlate, fSelectedSector));
625   }
626 }
627
628 //__________________________________________________________________
629 void AliTOFSDigitizer::SimulateDetectorResponse(Float_t z0, Float_t x0, Float_t geantTime, Int_t& nActivatedPads, Int_t& nFiredPads, Bool_t* isFired, Int_t* nPlace, Float_t* qInduced, Float_t* tofTime, Float_t& averageTime)
630 {
631   // Description:
632   // Input:  z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
633   //         geantTime - time generated by Geant, ns
634   // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
635   //         nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
636   //         qInduced[iPad]- charge induced on pad, arb. units
637   //                         this array is initialized at zero by the caller
638   //         tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
639   //                                   this array is initialized at zero by the caller
640   //         averageTime - time given by pad hited by the Geant track taking into account the times (weighted) given by the pads fired for edge effect also.
641   //                       The weight is given by the qInduced[iPad]/qCenterPad
642   //                                   this variable is initialized at zero by the caller
643   //         nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
644   //                                   this variable is initialized at zero by the caller
645   //
646   // Description of used variables:
647   //         eff[iPad] - efficiency of the pad
648   //         res[iPad] - resolution of the pad, ns
649   //         timeWalk[iPad] - time walk of the pad, ns
650   //         timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
651   //         PadId[iPad] - Pad Identifier
652   //                    E | F    -->   PadId[iPad] = 5 | 6
653   //                    A | B    -->   PadId[iPad] = 1 | 2
654   //                    C | D    -->   PadId[iPad] = 3 | 4
655   //         nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
656   //         qCenterPad - charge extimated for each pad, arb. units
657   //         weightsSum - sum of weights extimated for each pad fired, arb. units
658   
659   const Float_t kSigmaForTail[2] = {AliTOFGeometry::SigmaForTail1(),AliTOFGeometry::SigmaForTail2()}; //for tail                                                   
660   Int_t iz = 0, ix = 0;
661   Float_t dX = 0., dZ = 0., x = 0., z = 0.;
662   Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
663   Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
664   Float_t logOfqInd = 0.;
665   Float_t weightsSum = 0.;
666   Int_t nTail[4]  = {0,0,0,0};
667   Int_t padId[4]  = {0,0,0,0};
668   Float_t eff[4]  = {0.,0.,0.,0.};
669   Float_t res[4]  = {0.,0.,0.,0.};
670   //  Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
671   Float_t qCenterPad = 1.;
672   Float_t timeWalk[4]  = {0.,0.,0.,0.};
673   Float_t timeDelay[4] = {0.,0.,0.,0.};
674   
675   nActivatedPads = 0;
676   nFiredPads = 0;
677   
678   (z0 <= 0) ? iz = 0 : iz = 1;
679   dZ = z0 + (0.5 * AliTOFGeometry::NpadZ() - iz - 0.5) * AliTOFGeometry::ZPad(); // hit position in the pad frame, (0,0) - center of the pad
680   z = 0.5 * AliTOFGeometry::ZPad() - TMath::Abs(dZ);                               // variable for eff., res. and timeWalk. functions
681   iz++;                                                                              // z row: 1, ..., AliTOFGeometry::NpadZ = 2
682   ix = (Int_t)((x0 + 0.5 * AliTOFGeometry::NpadX() * AliTOFGeometry::XPad()) / AliTOFGeometry::XPad());
683   dX = x0 + (0.5 * AliTOFGeometry::NpadX() - ix - 0.5) * AliTOFGeometry::XPad(); // hit position in the pad frame, (0,0) - center of the pad
684   x = 0.5 * AliTOFGeometry::XPad() - TMath::Abs(dX);                               // variable for eff., res. and timeWalk. functions;
685   ix++;                                                                              // x row: 1, ..., AliTOFGeometry::NpadX = 48
686   
687   ////// Pad A:
688   nActivatedPads++;
689   nPlace[nActivatedPads-1] = (iz - 1) * AliTOFGeometry::NpadX() + ix;
690   qInduced[nActivatedPads-1] = qCenterPad;
691   padId[nActivatedPads-1] = 1;
692   
693   if (fEdgeEffect == 0) {
694     eff[nActivatedPads-1] = fEffCenter;
695     if (gRandom->Rndm() < eff[nActivatedPads-1]) {
696       nFiredPads = 1;
697       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + fResCenter * fResCenter); // ns
698       isFired[nActivatedPads-1] = kTRUE;
699       tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
700       averageTime = tofTime[nActivatedPads-1];
701     }
702   } else {
703      
704     if(z < h) {
705       if(z < h2) {
706         effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
707       } else {
708         effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
709       }
710       resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
711       timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
712       nTail[nActivatedPads-1] = 1;
713     } else {
714       effZ = fEffCenter;
715       resZ = fResCenter;
716       timeWalkZ = fTimeWalkCenter;
717     }
718     
719     if(x < h) {
720       if(x < h2) {
721         effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
722       } else {
723         effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
724       }
725       resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
726       timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
727       nTail[nActivatedPads-1] = 1;
728     } else {
729       effX = fEffCenter;
730       resX = fResCenter;
731       timeWalkX = fTimeWalkCenter;
732     }
733     
734     (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
735     (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns
736     (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 *  timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
737
738
739     ////// Pad B:
740     if(z < k2) {
741       effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
742     } else {
743       effZ = fEff3Boundary * (k - z) / (k - k2);
744     }
745     resZ = fResBoundary + fResSlope * z / k;
746     timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
747     
748     if(z < k && z > 0) {
749       if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
750         nActivatedPads++;
751         nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX();
752         eff[nActivatedPads-1] = effZ;
753         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns 
754         timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
755         nTail[nActivatedPads-1] = 2;
756         if (fTimeDelayFlag) {
757           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
758           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
759           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
760           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
761           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
762         } else {
763           timeDelay[nActivatedPads-1] = 0.;
764         }
765         padId[nActivatedPads-1] = 2;
766       }
767     }
768
769     
770     ////// Pad C, D, E, F:
771     if(x < k2) {
772       effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
773     } else {
774       effX = fEff3Boundary * (k - x) / (k - k2);
775     }
776     resX = fResBoundary + fResSlope*x/k;
777     timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
778     
779     if(x < k && x > 0) {
780       //   C:
781       if(ix > 1 && dX < 0) {
782         nActivatedPads++;
783         nPlace[nActivatedPads-1] = nPlace[0] - 1;
784         eff[nActivatedPads-1] = effX;
785         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX); // ns 
786         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
787         nTail[nActivatedPads-1] = 2;
788         if (fTimeDelayFlag) {
789           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
790           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
791           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
792           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
793           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
794         } else {
795           timeDelay[nActivatedPads-1] = 0.;
796         }
797         padId[nActivatedPads-1] = 3;
798
799         //     D:
800         if(z < k && z > 0) {
801           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
802             nActivatedPads++;
803             nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() - 1;
804             eff[nActivatedPads-1] = effX * effZ;
805             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns
806             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
807             
808             nTail[nActivatedPads-1] = 2;
809             if (fTimeDelayFlag) {
810               if (TMath::Abs(x) < TMath::Abs(z)) {
811                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
812                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
813                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
814                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
815               } else {
816                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
817                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
818                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
819                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
820               }
821               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
822             } else {
823               timeDelay[nActivatedPads-1] = 0.;
824             }
825             padId[nActivatedPads-1] = 4;
826           }
827         }  // end D
828       }  // end C
829       
830       //   E:
831       if(ix < AliTOFGeometry::NpadX() && dX > 0) {
832         nActivatedPads++;
833         nPlace[nActivatedPads-1] = nPlace[0] + 1;
834         eff[nActivatedPads-1] = effX;
835         res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(fAddTRes*fAddTRes + resX * resX)); // ns
836         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
837         nTail[nActivatedPads-1] = 2;
838         if (fTimeDelayFlag) {
839           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
840           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
841           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
842           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
843           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
844         } else {
845           timeDelay[nActivatedPads-1] = 0.;
846         }
847         padId[nActivatedPads-1] = 5;
848
849
850         //     F:
851         if(z < k && z > 0) {
852           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
853             nActivatedPads++;
854             nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() + 1;
855             eff[nActivatedPads - 1] = effX * effZ;
856             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns
857             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
858             nTail[nActivatedPads-1] = 2;
859             if (fTimeDelayFlag) {
860               if (TMath::Abs(x) < TMath::Abs(z)) {
861                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
862                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
863                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
864                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
865               } else {
866                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
867                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
868                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
869                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
870               }
871               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
872             } else {
873               timeDelay[nActivatedPads-1] = 0.;
874             }
875             padId[nActivatedPads-1] = 6;
876           }
877         }  // end F
878       }  // end E
879     } // end if(x < k)
880
881
882     for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
883       if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
884       if(gRandom->Rndm() < eff[iPad]) {
885         isFired[iPad] = kTRUE;
886         nFiredPads++;
887         if(fEdgeTails) {
888           if(nTail[iPad] == 0) {
889             tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
890           } else {
891             ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
892             Double_t timeAB = ftail->GetRandom();
893             tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
894           }
895         } else {
896           tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
897         }
898         if (fAverageTimeFlag) {
899           averageTime += tofTime[iPad] * qInduced[iPad];
900           weightsSum += qInduced[iPad];
901         } else {
902           averageTime += tofTime[iPad];
903           weightsSum += 1.;
904         }
905       }
906     }
907     if (weightsSum!=0) averageTime /= weightsSum;
908   } // end else (fEdgeEffect != 0)
909 }
910
911 //__________________________________________________________________
912 void AliTOFSDigitizer::PrintParameters()const
913 {
914   //
915   // Print parameters used for sdigitization
916   //
917   AliInfo(Form(" ------------------- %s -------------", GetName()));
918   AliInfo(" Parameters used for TOF SDigitization ");
919   //  Printing the parameters
920   
921   AliInfo(Form(" Number of events:                       %i ", (fEvent2-fEvent1)));
922   AliInfo(Form(" from event %i to event %i", fEvent1, (fEvent2-1)));
923   AliInfo(Form(" Time Resolution (ns) %d  Pad Efficiency: %d ", fTimeResolution, fpadefficiency));
924   AliInfo(Form(" Edge Effect option:  %d", fEdgeEffect));
925
926   AliInfo(" Boundary Effect Simulation Parameters ");
927   AliInfo(Form(" Hparameter: %d  H2parameter: %d  Kparameter: %d  K2parameter: %d", fHparameter, fH2parameter, fKparameter, fK2parameter));
928   AliInfo(Form(" Efficiency in the central region of the pad: %d", fEffCenter));
929   AliInfo(Form(" Efficiency at the boundary region of the pad: %d", fEffBoundary));
930   AliInfo(Form(" Efficiency value at H2parameter %d", fEff2Boundary));
931   AliInfo(Form(" Efficiency value at K2parameter %d", fEff3Boundary));
932   AliInfo(Form(" Resolution (ps) in the central region of the pad: %d", fResCenter));
933   AliInfo(Form(" Resolution (ps) at the boundary of the pad      : %d", fResBoundary));
934   AliInfo(Form(" Slope (ps/K) for neighbouring pad               : %d", fResSlope));
935   AliInfo(Form(" Time walk (ps) in the central region of the pad : %d", fTimeWalkCenter));
936   AliInfo(Form(" Time walk (ps) at the boundary of the pad       : %d", fTimeWalkBoundary));
937   AliInfo(Form(" Slope (ps/K) for neighbouring pad               : %d", fTimeWalkSlope));
938   AliInfo(" Pulse Heigth Simulation Parameters ");
939   AliInfo(Form(" Flag for delay due to the PulseHeightEffect  : %d", fTimeDelayFlag));
940   AliInfo(Form(" Pulse Height Slope                           : %d", fPulseHeightSlope));
941   AliInfo(Form(" Time Delay Slope                             : %d", fTimeDelaySlope));
942   AliInfo(Form(" Minimum charge amount which could be induced : %d", fMinimumCharge));
943   AliInfo(Form(" Smearing in charge in (q1/q2) vs x plot      : %d", fChargeSmearing));
944   AliInfo(Form(" Smearing in log of charge ratio              : %d", fLogChargeSmearing));
945   AliInfo(Form(" Smearing in time in time vs log(q1/q2) plot  : %d", fTimeSmearing));
946   AliInfo(Form(" Flag for average time                        : %d", fAverageTimeFlag));
947   AliInfo(Form(" Edge tails option                            : %d", fEdgeTails));
948   
949 }