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