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