- // cout << "hit: " << ii+1 << " of " << nhits << endl;
- AliITShit *hit = (AliITShit*) fHits->At(ii);
- AliITShit *hit1 = 0;
-
- // Take into account all hits when several GEANT steps are carried out
- // inside the silicon
- // Get and use the status of hit(track):
- // 66 - for entering hit,
- // 65 - for inside hit,
- // 68 - for exiting hit,
- // 33 - for stopping hit.
-
- //Int_t status = hit->GetTrackStatus();
- Int_t status1 = 0;
- Int_t hitDetector = hit->GetDetector();
- Float_t depEnergy = 0.;
- if(hit->StatusEntering()) { // to be coupled to following hit
- idhit=ii;
- hit->GetPositionL(xL[0],xL[1],xL[2]);
- if(ii<nhits-1) ii++;
- hit1 = (AliITShit*) fHits->At(ii);
- hit1->GetPositionL(xL1[0],xL1[1],xL1[2]);
- status1 = hit1->GetTrackStatus();
- depEnergy = kconv*hit1->GetIonization();
- } else {
- depEnergy = kconv*hit->GetIonization(); // Deposited energy in keV
- hit->GetPositionL(xL1[0],xL1[1],xL1[2]);
- }
- // cout << "status: " << status << ", status1: " << status1 << ", dE: " << depEnergy << endl;
- if(fFlag && status1 == 33) continue;
-
- Int_t nOfSplits = 1;
-
- // hit->Print();
-
-// Int_t status1 = -1;
-// Int_t ctr = 0;
- //Take now the entering and inside hits only
-// if(status == 66) {
-// do {
-// if(ii<nhits-1) ii++;
-// hit1 = (AliITShit*) fHits->At(ii);
-// hit1->GetPositionL(xL1[0],xL1[1],xL1[2]);
-// status1 = hit1->GetTrackStatus();
-// depEnergy += kconv*hit1->GetIonization();
-// if(fFlag && status1 == 65) ctr++;
-// } while(status1 != 68 && status1 != 33);
-// }
-
-
- // scale path to simulate a perpendicular track
- // continue if the particle did not lose energy
- // passing through detector
- if (!depEnergy) {
- printf("This particle has passed without losing energy!\n");
- continue;
- }
- Float_t pathInSDD = TMath::Sqrt((xL[0]-xL1[0])*(xL[0]-xL1[0])+(xL[1]-xL1[1])*(xL[1]-xL1[1])+(xL[2]-xL1[2])*(xL[2]-xL1[2]));
-
- if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
- Float_t Drft = (xL1[0]+xL[0])*0.5;
- Float_t drPath = 10000.*Drft;
- if(drPath < 0) drPath = -drPath;
- drPath = sddLength-drPath;
- if(drPath < 0) {
- cout << "Warning: negative drift path " << drPath << endl;
- continue;
- }
-
- // Drift Time
- Float_t drTime = drPath/driftSpeed;
- // Signal 2d Shape
- Float_t dfCoeff, s1;
- fResponse->DiffCoeff(dfCoeff,s1);
-
- // Squared Sigma along the anodes
- Double_t sig2A = 2.*dfCoeff*drTime+s1*s1;
- Double_t sigA = TMath::Sqrt(sig2A);
- if(pathInSDD) {
- nOfSplits = (Int_t) (1 + 10000.*pathInSDD/sigA);
- //cout << "nOfSplits: " << nOfSplits << ", sigA: " << sigA << ", path: " << pathInSDD << endl;
- }
- if(fFlag) nOfSplits = 1;
- depEnergy /= nOfSplits;
-
- for(Int_t kk=0;kk<nOfSplits;kk++) {
- Float_t avDrft =
- xL[0]+(xL1[0]-xL[0])*((kk+0.5)/((Float_t) nOfSplits));
- Float_t avAnode =
- xL[2]+(xL1[2]-xL[2])*((kk+0.5)/((Float_t) nOfSplits));
- Float_t driftPath = 10000.*avDrft;
-
- Int_t iWing = 2;
- if(driftPath < 0) {
- iWing = 1;
- driftPath = -driftPath;
- }
- driftPath = sddLength-driftPath;
- Int_t detector = 2*(hitDetector-1) + iWing;
- if(driftPath < 0) {
- cout << "Warning: negative drift path " << driftPath << endl;
- continue;
- }
-
- // Drift Time
- Float_t driftTime = driftPath/driftSpeed;
- Int_t timeSample = (Int_t) (fScaleSize*driftTime/timeStep + 1);
- if(timeSample > fScaleSize*fMaxNofSamples) {
- cout << "Warning: Wrong Time Sample: " << timeSample << endl;
- continue;
- }
-
- // Anode
- Float_t xAnode = 10000.*(avAnode)/anodePitch + nofAnodes/2; // +1?
- if(xAnode*anodePitch > sddWidth || xAnode*anodePitch < 0.)
- { cout << "Warning: Z = " << xAnode*anodePitch << endl; }
- Int_t iAnode = (Int_t) (1.+xAnode); // xAnode?
- if(iAnode < 1 || iAnode > nofAnodes) {
- cout << "Warning: Wrong iAnode: " << iAnode << endl;
- continue;
- }
-
- // work with the idtrack=entry number in the TreeH for the moment
- //Int_t idhit,idtrack;
- //mod->GetHitTrackAndHitIndex(ii,idtrack,idhit);
- //Int_t idtrack=mod->GetHitTrackIndex(ii);
- // or store straight away the particle position in the array
- // of particles and take idhit=ii only when part is entering (this
- // requires FillModules() in the macro for analysis) :
- Int_t itrack = hit->GetTrack();
-
- // Signal 2d Shape
- Float_t diffCoeff, s0;
- fResponse->DiffCoeff(diffCoeff,s0);
-
- // Squared Sigma along the anodes
- Double_t sigma2A = 2.*diffCoeff*driftTime+s0*s0;
- Double_t sigmaA = TMath::Sqrt(sigma2A);
- Double_t sigmaT = sigmaA/driftSpeed;
- // Peak amplitude in nanoAmpere
- Double_t eVpairs = 3.6;
- Double_t amplitude = fScaleSize*160.*depEnergy/(timeStep*eVpairs*2.*acos(-1.)*sigmaT*sigmaA);
- amplitude *= timeStep/25.; // WARNING!!!!! Amplitude scaling to account for clock variations (reference value: 40 MHz)
- Double_t chargeloss = 1.-CHloss*driftPath/1000;
- amplitude *= chargeloss;
- Float_t nsigma=fResponse->NSigmaIntegration();
- Int_t nlookups = fResponse->GausNLookUp();
- Float_t width = 2.*nsigma/(nlookups-1);
- // Spread the charge
- // Pixel index
- Int_t ja = iAnode;
- Int_t jt = timeSample;
- Int_t ndiv = 2;
- Float_t nmul = 3.;
- if(driftTime > 1200.) {
- ndiv = 4;
- nmul = 1.5;
- }
- // Sub-pixel index
- Int_t nsplit = 4; // hard-wired
- nsplit = (nsplit+1)/2*2;
- // Sub-pixel size
- Double_t aStep = anodePitch/(nsplit*fScaleSize);
- Double_t tStep = timeStep/(nsplit*fScaleSize);
- // Define SDD window corresponding to the hit
- Int_t anodeWindow = (Int_t) (fScaleSize*nsigma*sigmaA/anodePitch + 1);
- Int_t timeWindow = (Int_t) (fScaleSize*nsigma*sigmaT/timeStep + 1);
- Int_t jamin = (ja - anodeWindow/ndiv - 1)*fScaleSize*nsplit + 1;
- Int_t jamax = (ja + anodeWindow/ndiv)*fScaleSize*nsplit;
- if(jamin <= 0) jamin = 1;
- if(jamax > fScaleSize*nofAnodes*nsplit) jamax = fScaleSize*nofAnodes*nsplit;
- Int_t jtmin = (Int_t) (jt - timeWindow*nmul - 1)*nsplit + 1; //hard-wired
- Int_t jtmax = (Int_t) (jt + timeWindow*nmul)*nsplit; //hard-wired
- if(jtmin <= 0) jtmin = 1;
- if(jtmax > fScaleSize*fMaxNofSamples*nsplit) jtmax = fScaleSize*fMaxNofSamples*nsplit;
-
- // Spread the charge in the anode-time window
- Int_t ka;
- //cout << "jamin: " << jamin << ", jamax: " << jamax << endl;
- //cout << "jtmin: " << jtmin << ", jtmax: " << jtmax << endl;
- for(ka=jamin; ka <=jamax; ka++) {
- Int_t ia = (ka-1)/(fScaleSize*nsplit) + 1;
- if(ia <= 0) { cout << "Warning: ia < 1: " << endl; continue; }
- if(ia > nofAnodes) ia = nofAnodes;
- Double_t aExpo = (aStep*(ka-0.5)-xAnode*anodePitch)/sigmaA;
- Double_t anodeAmplitude = 0;
- if(TMath::Abs(aExpo) > nsigma) {
- anodeAmplitude = 0.;
- //cout << "aExpo: " << aExpo << endl;
- } else {
- Int_t i = (Int_t) ((aExpo+nsigma)/width);
- //cout << "eval ampl: " << i << ", " << amplitude << endl;
- anodeAmplitude = amplitude*fResponse->GausLookUp(i);
- //cout << "ampl: " << anodeAmplitude << endl;