Updated AliTOFSDigitizer
authorvicinanz <vicinanz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Apr 2002 10:06:34 +0000 (10:06 +0000)
committervicinanz <vicinanz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Apr 2002 10:06:34 +0000 (10:06 +0000)
TOF/AliTOFSDigitizer.cxx
TOF/AliTOFSDigitizer.h

index a87eef6..6af0712 100644 (file)
@@ -29,6 +29,7 @@
 
 #include "AliTOFHitMap.h"
 #include "AliTOFSDigit.h"
+#include "AliTOFConstants.h"
 #include "AliTOFhit.h"
 #include "AliTOF.h"
 #include "AliTOFv1.h"
@@ -46,6 +47,7 @@
 #include "TSystem.h"
 #include "TROOT.h"
 #include "TFolder.h"
+#include <TF1.h>
 #include <stdlib.h>
 #include <iostream.h>
 #include <fstream.h>
@@ -59,13 +61,18 @@ ClassImp(AliTOFSDigitizer)
   fNevents = 0 ;     
 //  fSDigits = 0 ;
   fHits = 0 ;
-
+  ftail    = 0;
 }
            
 //____________________________________________________________________________ 
   AliTOFSDigitizer::AliTOFSDigitizer(char* HeaderFile,char *SdigitsFile ):TTask("AliTOFSDigitizer","") 
 {
   fNevents = 0 ;     // Number of events to digitize, 0 means all evens in current file
+  ftail    = 0;
+
+  // init parameters for sdigitization
+  InitParameters();
+
   // add Task to //root/Tasks folder
   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
   roottasks->Add(this) ; 
@@ -75,8 +82,65 @@ ClassImp(AliTOFSDigitizer)
   AliTOFSDigitizer::~AliTOFSDigitizer()
 {
   // dtor
+  if (ftail)
+    {
+      delete ftail;
+      ftail = 0;
+    }
+}
+
+//____________________________________________________________________________ 
+void AliTOFSDigitizer::InitParameters()
+{
+  // set parameters for detector simulation
+
+  fTimeResolution =0.120;
+  fpadefficiency  =0.99 ;
+  fEdgeEffect     = 2   ;
+  fEdgeTails      = 0   ;
+  fHparameter     = 0.4 ;
+  fH2parameter    = 0.15;
+  fKparameter     = 0.5 ;
+  fK2parameter    = 0.35;
+  fEffCenter      = fpadefficiency;
+  fEffBoundary    = 0.65;
+  fEff2Boundary   = 0.90;
+  fEff3Boundary   = 0.08;
+  fResCenter      = 50. ;
+  fResBoundary    = 70. ;
+  fResSlope       = 40. ;
+  fTimeWalkCenter = 0.  ;
+  fTimeWalkBoundary=0.  ;
+  fTimeWalkSlope  = 0.  ;
+  fTimeDelayFlag  = 1   ;
+  fPulseHeightSlope=2.0 ;
+  fTimeDelaySlope =0.060;
+  // was fMinimumCharge = TMath::Exp(fPulseHeightSlope*fKparameter/2.);
+  fMinimumCharge = TMath::Exp(-fPulseHeightSlope*fHparameter);
+  fChargeSmearing=0.0   ;
+  fLogChargeSmearing=0.13;
+  fTimeSmearing   =0.022;
+  fAverageTimeFlag=0    ;
+
+}
+
+//__________________________________________________________________
+Double_t TimeWithTail(Double_t* x, Double_t* par)
+{
+  // sigma - par[0], alpha - par[1], part - par[2]
+  //  at x<part*sigma - gauss
+  //  at x>part*sigma - TMath::Exp(-x/alpha)
+  Float_t xx =x[0];
+  Double_t f;
+  if(xx<par[0]*par[2]) {
+    f = TMath::Exp(-xx*xx/(2*par[0]*par[0]));
+  } else {
+    f = TMath::Exp(-(xx-par[0]*par[2])/par[1]-0.5*par[2]*par[2]);
+  }
+  return f;
 }
 
+
 //____________________________________________________________________________
 void AliTOFSDigitizer::Exec(Option_t *option) { 
 
@@ -88,6 +152,8 @@ void AliTOFSDigitizer::Exec(Option_t *option) {
     return;
   }
 
+  if (fEdgeTails) ftail = new TF1("tail",TimeWithTail,-2,2,3);
+
   if (fNevents == 0)
     fNevents = (Int_t) gAlice->TreeE()->GetEntries();
 
@@ -114,7 +180,7 @@ void AliTOFSDigitizer::Exec(Option_t *option) {
     AliTOFhit *tofHit;
     TClonesArray *TOFhits = TOF->Hits();
 
-// create hit map
+    // create hit map
     AliTOFHitMap *hitMap = new AliTOFHitMap(TOF->SDigits());
 
     Int_t ntracks = static_cast<Int_t>(TH->GetEntries());
@@ -124,42 +190,64 @@ void AliTOFSDigitizer::Exec(Option_t *option) {
       TH->GetEvent(track);
       particle = gAlice->Particle(track);
       Int_t nhits = TOFhits->GetEntriesFast();
+      // cleaning all hits of the same track in the same pad volume
+      // it is a rare event, however it happens
+
+      Int_t previousTrack =0;
+      Int_t previousSector=0;
+      Int_t previousPlate =0;
+      Int_t previousStrip =0;
+      Int_t previousPadX  =0;
+      Int_t previousPadZ  =0;
 
       for (Int_t hit = 0; hit < nhits; hit++)
       {
        tofHit = (AliTOFhit *) TOFhits->UncheckedAt(hit);
+       Int_t tracknum = tofHit->GetTrack();
        vol[0] = tofHit->GetSector();
        vol[1] = tofHit->GetPlate();
        vol[2] = tofHit->GetStrip();
        vol[3] = tofHit->GetPadx();
        vol[4] = tofHit->GetPadz();
 
-       // 95% of efficiency to be inserted here
-       // edge effect to be inserted here
-       // cross talk  to be inserted here
-
-       Float_t idealtime = tofHit->GetTof(); // unit s
-       idealtime *= 1.E+12;  // conversion from s to ps
-       // fTimeRes is given usually in ps
-       Float_t tdctime   = gRandom->Gaus(idealtime, TOF->GetTimeRes());
-       digit[0] = tdctime;
-
-       // typical Landau Distribution to be inserted here
-       // instead of Gaussian Distribution
-       Float_t idealcharge = tofHit->GetEdep();
-       Float_t adccharge = gRandom->Gaus(idealcharge, TOF->GetChrgRes());
-       digit[1] = adccharge;
-       Int_t tracknum = tofHit->GetTrack();
+       Bool_t isCloneOfThePrevious=((tracknum==previousTrack) && (vol[0]==previousSector) && (vol[1]==previousPlate) && (vol[2]==previousStrip) && (vol[3]==previousPadX) && (vol[4]==previousPadZ));
+       
+       if(!isCloneOfThePrevious){
+         // update "previous" values
+         // in fact, we are yet in the future, so the present is past
+         previousTrack=tracknum;
+         previousSector=vol[0];
+         previousPlate=vol[1];
+         previousStrip=vol[2];
+         previousPadX=vol[3];
+         previousPadZ=vol[4];
 
-       // check if two digit are on the same pad; in that case we sum
-       // the two or more digits
-       if (hitMap->TestHit(vol) != kEmpty) {
-         AliTOFSDigit *sdig = static_cast<AliTOFSDigit*>(hitMap->GetHit(vol));
-         sdig->Update(tdctime,adccharge,tracknum);
-       } else {
-         TOF->AddSDigit(tracknum, vol, digit);
-         hitMap->SetHit(vol);
-       }
+         // 95% of efficiency to be inserted here
+         // edge effect to be inserted here
+         // cross talk  to be inserted here
+         
+         Float_t idealtime = tofHit->GetTof(); // unit s
+         idealtime *= 1.E+12;  // conversion from s to ps
+         // fTimeRes is given usually in ps
+         Float_t tdctime   = gRandom->Gaus(idealtime, TOF->GetTimeRes());
+         digit[0] = tdctime;
+         
+         // typical Landau Distribution to be inserted here
+         // instead of Gaussian Distribution
+         Float_t idealcharge = tofHit->GetEdep();
+         Float_t adccharge = gRandom->Gaus(idealcharge, TOF->GetChrgRes());
+         digit[1] = adccharge;
+         
+         // check if two digit are on the same pad; in that case we sum
+         // the two or more digits
+         if (hitMap->TestHit(vol) != kEmpty) {
+           AliTOFSDigit *sdig = static_cast<AliTOFSDigit*>(hitMap->GetHit(vol));
+           sdig->Update(tdctime,adccharge,tracknum);
+         } else {
+           TOF->AddSDigit(tracknum, vol, digit);
+           hitMap->SetHit(vol);
+         }
+       } // close if(!isCloneOfThePrevious)
       } // end loop on hits for the current track
     } // end loop on ntracks
 
@@ -189,3 +277,326 @@ void AliTOFSDigitizer::Print(Option_t* option)const
     cout << "    Writing SDigitis to file  " << (char*) fSDigitsFile.Data() << endl ;
 
 }
+
+//__________________________________________________________________
+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)
+{
+  // Description:
+  // Input:  z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
+  //         geantTime - time generated by Geant, ns
+  // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
+  //         nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
+  //         qInduced[iPad]- charge induced on pad, arb. units
+  //                         this array is initialized at zero by the caller
+  //         tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
+  //                                   this array is initialized at zero by the caller
+  //         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.
+  //                       The weight is given by the qInduced[iPad]/qCenterPad
+  //                                   this variable is initialized at zero by the caller
+  //         nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
+  //                                   this variable is initialized at zero by the caller
+  //
+  // Description of used variables:
+  //         eff[iPad] - efficiency of the pad
+  //         res[iPad] - resolution of the pad, ns
+  //         timeWalk[iPad] - time walk of the pad, ns
+  //         timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
+  //         PadId[iPad] - Pad Identifier
+  //                    E | F    -->   PadId[iPad] = 5 | 6
+  //                    A | B    -->   PadId[iPad] = 1 | 2
+  //                    C | D    -->   PadId[iPad] = 3 | 4
+  //         nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
+  //         qCenterPad - charge extimated for each pad, arb. units
+  //         weightsSum - sum of weights extimated for each pad fired, arb. units
+  
+  const Float_t kSigmaForTail[2] = {AliTOFConstants::fgkSigmaForTail1,AliTOFConstants::fgkSigmaForTail2}; //for tail                                                   
+  Int_t iz = 0, ix = 0;
+  Float_t dX = 0., dZ = 0., x = 0., z = 0.;
+  Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
+  Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
+  Float_t logOfqInd = 0.;
+  Float_t weightsSum = 0.;
+  Int_t nTail[4]  = {0,0,0,0};
+  Int_t padId[4]  = {0,0,0,0};
+  Float_t eff[4]  = {0.,0.,0.,0.};
+  Float_t res[4]  = {0.,0.,0.,0.};
+  //  Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
+  Float_t qCenterPad = 1.;
+  Float_t timeWalk[4]  = {0.,0.,0.,0.};
+  Float_t timeDelay[4] = {0.,0.,0.,0.};
+  
+  nActivatedPads = 0;
+  nFiredPads = 0;
+  
+  (z0 <= 0) ? iz = 0 : iz = 1;
+  dZ = z0 + (0.5 * AliTOFConstants::fgkNpadZ - iz - 0.5) * AliTOFConstants::fgkZPad; // hit position in the pad frame, (0,0) - center of the pad
+  z = 0.5 * AliTOFConstants::fgkZPad - TMath::Abs(dZ);                               // variable for eff., res. and timeWalk. functions
+  iz++;                                                                              // z row: 1, ..., AliTOFConstants::fgkNpadZ = 2
+  ix = (Int_t)((x0 + 0.5 * AliTOFConstants::fgkNpadX * AliTOFConstants::fgkXPad) / AliTOFConstants::fgkXPad);
+  dX = x0 + (0.5 * AliTOFConstants::fgkNpadX - ix - 0.5) * AliTOFConstants::fgkXPad; // hit position in the pad frame, (0,0) - center of the pad
+  x = 0.5 * AliTOFConstants::fgkXPad - TMath::Abs(dX);                               // variable for eff., res. and timeWalk. functions;
+  ix++;                                                                              // x row: 1, ..., AliTOFConstants::fgkNpadX = 48
+  
+  ////// Pad A:
+  nActivatedPads++;
+  nPlace[nActivatedPads-1] = (iz - 1) * AliTOFConstants::fgkNpadX + ix;
+  qInduced[nActivatedPads-1] = qCenterPad;
+  padId[nActivatedPads-1] = 1;
+  
+  if (fEdgeEffect == 0) {
+    eff[nActivatedPads-1] = fEffCenter;
+    if (gRandom->Rndm() < eff[nActivatedPads-1]) {
+      nFiredPads = 1;
+      res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns;
+      isFired[nActivatedPads-1] = kTRUE;
+      tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
+      averageTime = tofTime[nActivatedPads-1];
+    }
+  } else {
+     
+    if(z < h) {
+      if(z < h2) {
+       effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
+      } else {
+       effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
+      }
+      resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
+      timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
+      nTail[nActivatedPads-1] = 1;
+    } else {
+      effZ = fEffCenter;
+      resZ = fResCenter;
+      timeWalkZ = fTimeWalkCenter;
+    }
+    
+    if(x < h) {
+      if(x < h2) {
+       effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
+      } else {
+       effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
+      }
+      resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
+      timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
+      nTail[nActivatedPads-1] = 1;
+    } else {
+      effX = fEffCenter;
+      resX = fResCenter;
+      timeWalkX = fTimeWalkCenter;
+    }
+    
+    (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
+    (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns
+    (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 *  timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
+
+
+    ////// Pad B:
+    if(z < k2) {
+      effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
+    } else {
+      effZ = fEff3Boundary * (k - z) / (k - k2);
+    }
+    resZ = fResBoundary + fResSlope * z / k;
+    timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
+    
+    if(z < k && z > 0) {
+      if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
+       nActivatedPads++;
+       nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX;
+       eff[nActivatedPads-1] = effZ;
+       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns 
+       timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
+       nTail[nActivatedPads-1] = 2;
+       if (fTimeDelayFlag) {
+         //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
+         //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
+         qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
+         logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
+         timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
+       } else {
+         timeDelay[nActivatedPads-1] = 0.;
+       }
+       padId[nActivatedPads-1] = 2;
+      }
+    }
+
+    
+    ////// Pad C, D, E, F:
+    if(x < k2) {
+      effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
+    } else {
+      effX = fEff3Boundary * (k - x) / (k - k2);
+    }
+    resX = fResBoundary + fResSlope*x/k;
+    timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
+    
+    if(x < k && x > 0) {
+      //   C:
+      if(ix > 1 && dX < 0) {
+       nActivatedPads++;
+       nPlace[nActivatedPads-1] = nPlace[0] - 1;
+       eff[nActivatedPads-1] = effX;
+       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns 
+       timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
+       nTail[nActivatedPads-1] = 2;
+       if (fTimeDelayFlag) {
+         //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
+         //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
+         qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
+         logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
+         timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
+       } else {
+         timeDelay[nActivatedPads-1] = 0.;
+       }
+       padId[nActivatedPads-1] = 3;
+
+       //     D:
+       if(z < k && z > 0) {
+         if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
+           nActivatedPads++;
+           nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX - 1;
+           eff[nActivatedPads-1] = effX * effZ;
+           (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
+           (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
+           
+           nTail[nActivatedPads-1] = 2;
+           if (fTimeDelayFlag) {
+             if (TMath::Abs(x) < TMath::Abs(z)) {
+               //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
+               //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
+               qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
+               logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
+             } else {
+               //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
+               //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
+               qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
+               logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
+             }
+             timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
+           } else {
+             timeDelay[nActivatedPads-1] = 0.;
+           }
+           padId[nActivatedPads-1] = 4;
+         }
+       }  // end D
+      }  // end C
+      
+      //   E:
+      if(ix < AliTOFConstants::fgkNpadX && dX > 0) {
+       nActivatedPads++;
+       nPlace[nActivatedPads-1] = nPlace[0] + 1;
+       eff[nActivatedPads-1] = effX;
+       res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns
+       timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
+       nTail[nActivatedPads-1] = 2;
+       if (fTimeDelayFlag) {
+         //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
+         //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
+         qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
+         logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
+         timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
+       } else {
+         timeDelay[nActivatedPads-1] = 0.;
+       }
+       padId[nActivatedPads-1] = 5;
+
+
+       //     F:
+       if(z < k && z > 0) {
+         if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
+           nActivatedPads++;
+           nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX + 1;
+           eff[nActivatedPads - 1] = effX * effZ;
+           (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
+           (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
+           nTail[nActivatedPads-1] = 2;
+           if (fTimeDelayFlag) {
+             if (TMath::Abs(x) < TMath::Abs(z)) {
+               //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
+               //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
+               qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
+               logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
+             } else {
+               //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
+               //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
+               qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
+               logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
+             }
+             timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
+           } else {
+             timeDelay[nActivatedPads-1] = 0.;
+           }
+           padId[nActivatedPads-1] = 6;
+         }
+       }  // end F
+      }  // end E
+    } // end if(x < k)
+
+
+    for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
+      if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
+      if(gRandom->Rndm() < eff[iPad]) {
+       isFired[iPad] = kTRUE;
+       nFiredPads++;
+       if(fEdgeTails) {
+         if(nTail[iPad] == 0) {
+           tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
+         } else {
+           ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
+           Double_t timeAB = ftail->GetRandom();
+           tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
+         }
+       } else {
+         tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
+       }
+       if (fAverageTimeFlag) {
+         averageTime += tofTime[iPad] * qInduced[iPad];
+         weightsSum += qInduced[iPad];
+       } else {
+         averageTime += tofTime[iPad];
+         weightsSum += 1.;
+       }
+      }
+    }
+    if (weightsSum!=0) averageTime /= weightsSum;
+  } // end else (fEdgeEffect != 0)
+}
+
+//__________________________________________________________________
+void AliTOFSDigitizer::PrintParameters()const
+{
+  //
+  // Print parameters used for sdigitization
+  //
+  cout << " ------------------- "<< GetName() << " -------------" << endl ;
+  cout << " Parameters used for TOF SDigitization " << endl ;
+  //  Printing the parameters
+  
+  cout << " Number of events:                        " << fNevents << endl; 
+
+  cout << " Time Resolution (ns) "<< fTimeResolution <<" Pad Efficiency: "<< fpadefficiency << endl;
+  cout << " Edge Effect option:  "<<  fEdgeEffect<< endl;
+
+  cout << " Boundary Effect Simulation Parameters " << endl;
+  cout << " Hparameter: "<< fHparameter<<"  H2parameter:"<< fH2parameter <<"  Kparameter:"<< fKparameter<<"  K2parameter: "<< fK2parameter << endl;
+  cout << " Efficiency in the central region of the pad: "<< fEffCenter << endl;
+  cout << " Efficiency at the boundary region of the pad: "<< fEffBoundary << endl;
+  cout << " Efficiency value at H2parameter "<< fEff2Boundary << endl;
+  cout << " Efficiency value at K2parameter "<< fEff3Boundary << endl;
+  cout << " Resolution (ps) in the central region of the pad: "<< fResCenter << endl;
+  cout << " Resolution (ps) at the boundary of the pad      : "<< fResBoundary << endl;
+  cout << " Slope (ps/K) for neighbouring pad               : "<< fResSlope <<endl;
+  cout << " Time walk (ps) in the central region of the pad : "<< fTimeWalkCenter << endl;
+  cout << " Time walk (ps) at the boundary of the pad       : "<< fTimeWalkBoundary<< endl;
+  cout << " Slope (ps/K) for neighbouring pad               : "<< fTimeWalkSlope<<endl;
+  cout << " Pulse Heigth Simulation Parameters " << endl;
+  cout << " Flag for delay due to the PulseHeightEffect: "<< fTimeDelayFlag <<endl;
+  cout << " Pulse Height Slope                           : "<< fPulseHeightSlope<<endl;
+  cout << " Time Delay Slope                             : "<< fTimeDelaySlope<<endl;
+  cout << " Minimum charge amount which could be induced : "<< fMinimumCharge<<endl;
+  cout << " Smearing in charge in (q1/q2) vs x plot      : "<< fChargeSmearing<<endl;
+  cout << " Smearing in log of charge ratio              : "<< fLogChargeSmearing<<endl;
+  cout << " Smearing in time in time vs log(q1/q2) plot  : "<< fTimeSmearing<<endl;
+  cout << " Flag for average time                        : "<< fAverageTimeFlag<<endl;
+  cout << " Edge tails option                            : "<< fEdgeTails << endl;
+  
+}
index 9ec6e56..c871f74 100644 (file)
@@ -15,6 +15,8 @@
 #include "AliTOF.h"
 #include "AliDetector.h"
 
+class TF1;
+
 class AliTOFSDigitizer: public TTask {
 
 public:
@@ -22,7 +24,6 @@ public:
   AliTOFSDigitizer(char* HeaderFile, char *SdigitsFile = 0) ; 
 
   virtual ~AliTOFSDigitizer() ; // dtor
-  // Int_t    Digitize(Float_t Energy);
 
 //  char *GetSDigitsFile() const {return const_cast<char*>(fSDigitsFile.Data());}  
   const char *GetSDigitsFile() const {return fSDigitsFile.Data();}  
@@ -30,19 +31,114 @@ public:
   void  SetNEvents(Int_t Nevents) {fNevents = Nevents;}
   Int_t GetNEvents() const {return fNevents;}
   void SetSDigitsFile(char * file ) ;
+  void InitParameters();
+  virtual void PrintParameters() const ;
+  virtual void  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);
   virtual void Print(Option_t* option) const ;
   TClonesArray *SDigits() const {return fSDigits;}
   TClonesArray *Hits() const {return fHits;}
 
-
+  // setters and getters for detector simulation
+  // it summarizes all it is known about TOF strip 
+  void  SetPadefficiency(Float_t padefficiency)      {fpadefficiency=padefficiency;}
+  void  SetEdgeEffect(Int_t   edgeEffect)            {fEdgeEffect=edgeEffect;}
+  void  SetEdgeTails(Int_t   edgeTails)              {fEdgeTails=edgeTails;}
+  void  SetHparameter(Float_t hparameter)            {fHparameter=hparameter;}
+  void  SetH2parameter(Float_t h2parameter)          {fH2parameter=h2parameter;}
+  void  SetKparameter(Float_t kparameter)            {fKparameter=kparameter;}
+  void  SetK2parameter(Float_t k2parameter)          {fK2parameter=k2parameter;}
+  void  SetEffCenter(Float_t effCenter)              {fEffCenter=effCenter;}
+  void  SetEffBoundary(Float_t effBoundary)          {fEffBoundary=effBoundary;}
+  void  SetEff2Boundary(Float_t eff2Boundary)        {fEff2Boundary=eff2Boundary;}
+  void  SetEff3Boundary(Float_t eff3Boundary)        {fEff3Boundary=eff3Boundary;}
+  void  SetResCenter (Float_t resCenter)             {fResCenter=resCenter;}
+  void  SetResBoundary(Float_t resBoundary)          {fResBoundary=resBoundary;}
+  void  SetResSlope(Float_t resSlope)                {fResSlope=resSlope;}
+  void  SetTimeWalkCenter(Float_t timeWalkCenter)    {fTimeWalkCenter=timeWalkCenter;}
+  void  SetTimeWalkBoundary(Float_t timeWalkBoundary){fTimeWalkBoundary=timeWalkBoundary;}
+  void  SetTimeWalkSlope(Float_t timeWalkSlope)      {fTimeWalkSlope=timeWalkSlope;}
+
+  void  SetTimeDelayFlag(Int_t timeDelayFlag)        {fTimeDelayFlag=timeDelayFlag;}
+  void  SetPulseHeightSlope(Float_t pulseHeightSlope){fPulseHeightSlope=pulseHeightSlope;}
+  void  SetTimeDelaySlope(Float_t timeDelaySlope)    {fTimeDelaySlope=timeDelaySlope;}
+  void  SetMinimumCharge(Float_t minimumCharge)      {fMinimumCharge=minimumCharge;}
+  void  SetChargeSmearing(Float_t chargeSmearing)    {fChargeSmearing=chargeSmearing;}
+  void  SetLogChargeSmearing(Float_t logChargeSmearing){fLogChargeSmearing=logChargeSmearing;}
+  void  SetTimeSmearing(Float_t timeSmearing)        {fTimeSmearing=timeSmearing;}
+  void  SetAverageTimeFlag(Int_t averageTimeFlag)    {fAverageTimeFlag=averageTimeFlag;}
+
+  Float_t  GetPadefficiency()    const {return fpadefficiency;}
+  Int_t    GetEdgeEffect()       const {return fEdgeEffect;}
+  Int_t    GetEdgeTails()        const {return fEdgeTails;}
+  Float_t  GetHparameter()       const {return fHparameter;}
+  Float_t  GetH2parameter()      const {return fH2parameter;}
+  Float_t  GetKparameter()       const {return fKparameter;}
+  Float_t  GetK2parameter()      const {return fK2parameter;}
+  Float_t  GetEffCenter()        const {return fEffCenter;}
+  Float_t  GetEffBoundary()      const {return fEffBoundary;}
+  Float_t  GetEff2Boundary()     const {return fEff2Boundary;}
+  Float_t  GetEff3Boundary()     const {return fEff3Boundary;}
+  Float_t  GetResCenter ()       const {return fResCenter;}
+  Float_t  GetResBoundary()      const {return fResBoundary;}
+  Float_t  GetResSlope()         const {return fResSlope;}
+  Float_t  GetTimeWalkCenter()   const {return fTimeWalkCenter;}
+  Float_t  GetTimeWalkBoundary() const {return fTimeWalkBoundary;}
+  Float_t  GetTimeWalkSlope()    const {return fTimeWalkSlope;}
+  Int_t    GetTimeDelayFlag()    const {return fTimeDelayFlag;}
+  Float_t  GetPulseHeightSlope() const {return fPulseHeightSlope;}
+  Float_t  GetTimeDelaySlope()   const {return fTimeDelaySlope;}
+  Float_t  GetMinimumCharge()    const {return fMinimumCharge;}
+  Float_t  GetChargeSmearing()   const {return fChargeSmearing;}
+  Float_t  GetLogChargeSmearing()const {return fLogChargeSmearing;}
+  Float_t  GetTimeSmearing()     const {return fTimeSmearing;}
+  Int_t    GetAverageTimeFlag()  const {return fAverageTimeFlag;}
+  
 
 private:
   Int_t   fNevents;         // Number of events to digitize
   TString fSDigitsFile;     // output file 
   TClonesArray *fSDigits;   // array of summable digits
   TClonesArray *fHits;      // array of summable digits
+  TF1     *ftail;           // pointer to formula for time with tail
   TString fHeadersFile;     // input file
 
+  // detector response simulation
+  // Intrisic MRPC time resolution and pad (edge effect) parameters
+  Float_t fTimeResolution;  // time resolution of the MRPC (ns)
+  Float_t fpadefficiency;   // intrinsic pad efficiency, used if fEdgeEffect==0
+  Int_t   fEdgeEffect;      // edge effects option
+  Int_t   fEdgeTails;       // edge tails option
+  Float_t fHparameter;      // sensitive edge (to produce hits on the
+  // neighbouring pads) =0.7, new = 0.4 cm
+  Float_t fH2parameter;     // parameter to fit the efficiency
+  Float_t fKparameter;      // sensitive edge (going ahead towards the
+  // center no delay effects are suffered) =1.0, new = 0.5 cm
+  Float_t fK2parameter;     // parameter to fit the efficiency
+  // Pad Efficiency and Resolution parameters
+  Float_t fEffCenter;       // efficiency in the central region of the pad
+  Float_t fEffBoundary;     // efficiency at the boundary of the pad
+  Float_t fEff2Boundary;    // efficiency value at H2parameter
+  Float_t fEff3Boundary;    // efficiency value at K2parameter
+  Float_t fResCenter;       // resolution (ps) in the central region of the pad
+  Float_t fResBoundary;     // resolution (ps)  at the boundary of the pad
+  Float_t fResSlope;        // slope (ps/K) for neighbouring pad
+  // Time Walk parameters
+  Float_t fTimeWalkCenter;  // time walk (ps) in the central region of the pad
+  Float_t fTimeWalkBoundary;// time walk (ps) at the boundary of the pad
+  Float_t fTimeWalkSlope;   // slope (ps/K) for neighbouring pad
+  Int_t   fTimeDelayFlag;   // flag for delay due to the PulseHeightEffect
+  Float_t fPulseHeightSlope;// It determines the charge amount induced
+  // due to edge effect, using the formula
+  // qInduced=exp(-PulseHeightSlope*x)
+  Float_t fTimeDelaySlope;  // It determines the time delay. This is the slope
+  // in the T1-T2 vs log(q1/q2) plot
+  // ADC-TDC correlation parameters
+  Float_t fMinimumCharge;   // Minimum charge amount which could be induced
+  Float_t fChargeSmearing;  // Smearing in charge in (q1/q2) vs x plot
+  Float_t fLogChargeSmearing;// Smearing in log of charge ratio
+  Float_t fTimeSmearing;    // Smearing in time in time vs log(q1/q2) plot
+  Int_t   fAverageTimeFlag; // flag (see the setter for details)
+
  protected: