- Float_t nsigma=fResponse->NSigmaIntegration();
- // Spread the charge
- // Pixel index
- Int_t ja = iAnode;
- Int_t jt = timeSample;
- // 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/2 - 1)*fScaleSize*nsplit + 1;
- Int_t jamax = (ja + anodeWindow/2)*fScaleSize*nsplit;
- if(jamin <= 0) jamin = 1;
- if(jamax > fScaleSize*nofAnodes*nsplit) jamax = fScaleSize*nofAnodes*nsplit;
- Int_t jtmin = (jt - timeWindow*3 - 1)*nsplit + 1; //hard-wired
- Int_t jtmax = (jt + timeWindow*3)*nsplit; //hard-wired
- if(jtmin <= 0) jtmin = 1;
- if(jtmax > fScaleSize*fMaxNofSamples*nsplit) jtmax = fScaleSize*fMaxNofSamples*nsplit;
-
- Double_t rlAnode = log(aStep*amplitude);
-
- // Spread the charge in the anode-time window
- Int_t ka;
- 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 = rlAnode - 0.5*aExpo*aExpo;
- // Protect against overflows
- if(anodeAmplitude > -87.3)
- anodeAmplitude = exp(anodeAmplitude);
- else
- anodeAmplitude = 0;
- Int_t index = ((detector+1)%2)*nofAnodes+ia-1; // index starts from 0
- if(anodeAmplitude) {
- Double_t rlTime = log(tStep*anodeAmplitude);
- Int_t kt;
- for(kt=jtmin; kt<=jtmax; kt++) {
- Int_t it = (kt-1)/nsplit+1; // it starts from 1
- if(it<=0) { cout << "Warning: it < 1: " << endl; continue; }
- if(it>fScaleSize*fMaxNofSamples) it = fScaleSize*fMaxNofSamples;
- Double_t tExpo = (tStep*(kt-0.5)-driftTime)/sigmaT;
- Double_t timeAmplitude = rlTime - 0.5*tExpo*tExpo;
- // Protect against overflows
- if(timeAmplitude > -87.3){
- timeAmplitude = exp(timeAmplitude);
- } else
- timeAmplitude = 0;
-
- // build the list of digits for this module
- arg[0]=index;
- arg[1]=it;
- arg[2]=itrack;
- arg[3]=idhit;
- ListOfFiredCells(arg,timeAmplitude,list,padr);
- } // loop over time in window
- } // end if anodeAmplitude
- } // loop over anodes in window
- } // end loop over "sub-hits"
- } // end loop over hits
-
- // introduce the electronics effects and do zero-suppression if required
- Int_t nentries=list->GetEntriesFast();
- if (nentries) {
-
- //TStopwatch timer;
- ChargeToSignal();
- //timer.Stop(); timer.Print();
-
- const char *kopt=fResponse->ZeroSuppOption();
- ZeroSuppression(kopt);
- }
-
- // clean memory
- list->Delete();
- delete list;
-
- padr->Delete();
-
- fHitMap1->ClearMap();
- fHitMap2->ClearMap();
-
- //gObjectTable->Print();
-}
-
-
-//____________________________________________
-
-void AliITSsimulationSDD::ListOfFiredCells(Int_t *arg,Double_t timeAmplitude,
- TObjArray *list,TClonesArray *padr){
- // Returns the list of "fired" cells.
-
- Int_t index=arg[0];
- Int_t ik=arg[1];
- Int_t idtrack=arg[2];
- Int_t idhit=arg[3];
- Int_t counter=arg[4];
- Int_t countadr=arg[5];
-
- Double_t charge=timeAmplitude;
- charge += fHitMap2->GetSignal(index,ik-1);
- fHitMap2->SetHit(index, ik-1, charge);
-
- Int_t digits[3];
- Int_t it=(Int_t)((ik-1)/fScaleSize);
-
- digits[0]=index;
- digits[1]=it;
- digits[2]=(Int_t)timeAmplitude;
- Float_t phys;
- if (idtrack >= 0) phys=(Float_t)timeAmplitude;
- else phys=0;
-
- Double_t cellcharge=0.;
- AliITSTransientDigit* pdigit;
- // build the list of fired cells and update the info
- if (!fHitMap1->TestHit(index, it)) {
-
- new((*padr)[countadr++]) TVector(3);
- TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
- trinfo(0)=(Float_t)idtrack;
- trinfo(1)=(Float_t)idhit;
- trinfo(2)=(Float_t)timeAmplitude;
-
- list->AddAtAndExpand(
- new AliITSTransientDigit(phys,digits),counter);
-
- fHitMap1->SetHit(index, it, counter);
- counter++;
- pdigit=(AliITSTransientDigit*)list->
- At(list->GetLast());
- // list of tracks
- TObjArray *trlist=(TObjArray*)pdigit->TrackList();
- trlist->Add(&trinfo);
-
- } else {
- pdigit=
- (AliITSTransientDigit*) fHitMap1->GetHit(index, it);
- for(Int_t kk=0;kk<fScaleSize;kk++) {
- cellcharge += fHitMap2->GetSignal(index,fScaleSize*it+kk);
- }
- // update charge
- (*pdigit).fSignal=(Int_t)cellcharge;
- (*pdigit).fPhysics+=phys;
- // update list of tracks
- TObjArray* trlist=(TObjArray*)pdigit->TrackList();
- Int_t lastentry=trlist->GetLast();
- TVector *ptrkp=(TVector*)trlist->At(lastentry);
- TVector &trinfo=*ptrkp;
- Int_t lasttrack=Int_t(trinfo(0));
- //Int_t lasthit=Int_t(trinfo(1));
- Float_t lastcharge=(trinfo(2));
-
- if (lasttrack==idtrack ) {
- lastcharge+=(Float_t)timeAmplitude;
- trlist->RemoveAt(lastentry);
- trinfo(0)=lasttrack;
- //trinfo(1)=lasthit; // or idhit
- trinfo(1)=idhit;
- trinfo(2)=lastcharge;
- trlist->AddAt(&trinfo,lastentry);
- } else {
-
- new((*padr)[countadr++]) TVector(3);
- TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
- trinfo(0)=(Float_t)idtrack;
- trinfo(1)=(Float_t)idhit;
- trinfo(2)=(Float_t)timeAmplitude;
-
- trlist->Add(&trinfo);
- }
-
-#ifdef print
- // check the track list - debugging
- Int_t trk[20], htrk[20];
- Float_t chtrk[20];
- Int_t nptracks=trlist->GetEntriesFast();
- if (nptracks > 2) {
- Int_t tr;
- for (tr=0;tr<nptracks;tr++) {
- TVector *pptrkp=(TVector*)trlist->At(tr);
- TVector &pptrk=*pptrkp;
- trk[tr]=Int_t(pptrk(0));
- htrk[tr]=Int_t(pptrk(1));
- chtrk[tr]=(pptrk(2));
- printf("nptracks %d \n",nptracks);
- // set printings
- }
- } // end if nptracks
-#endif
- } // end if pdigit
-
- arg[4]=counter;
- arg[5]=countadr;
-
-
-}
-
-
-//____________________________________________
-
-void AliITSsimulationSDD::AddDigit(Int_t i, Int_t j, Int_t signal){
- // Adds a Digit.
- // tag with -1 signals coming from background tracks
- // tag with -2 signals coming from pure electronic noise
-
- Int_t digits[3], tracks[3], hits[3];
- Float_t phys, charges[3];
-
- Int_t trk[20], htrk[20];
- Float_t chtrk[20];
-
- Bool_t do10to8=fResponse->Do10to8();
-
- if(do10to8) signal=Convert8to10(signal);
- AliITSTransientDigit *obj = (AliITSTransientDigit*)fHitMap1->GetHit(i,j);
- digits[0]=i;
- digits[1]=j;
- digits[2]=signal;
- if (!obj) {
- phys=0;
- Int_t k;
- for (k=0;k<3;k++) {
- tracks[k]=-2;
- charges[k]=0;
- hits[k]=-1;
- }
- fITS->AddSimDigit(1,phys,digits,tracks,hits,charges);
- } else {
- phys=obj->fPhysics;
- TObjArray* trlist=(TObjArray*)obj->TrackList();
- Int_t nptracks=trlist->GetEntriesFast();
-
- if (nptracks > 20) {
- cout<<"Attention - nptracks > 20 "<<nptracks<<endl;
- nptracks=20;
+ // loop over path segments, init. some variables.
+ depEnergy /= nOfSplits;
+ nOfSplitsF = (Float_t) nOfSplits;
+ Float_t theAverage=0.,theSteps=0.;
+ for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
+ kkF = (Float_t) kk + 0.5;
+ avDrft = xL[0]+dxL[0]*kkF/nOfSplitsF;
+ avAnode = xL[2]+dxL[2]*kkF/nOfSplitsF;
+ theSteps+=1.;
+ theAverage+=avAnode;
+ zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
+ driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
+ driftSpeed+= fDetType->GetResponseSDD()->GetDeltaVDrift(fModule,zAnode>255);
+ driftPath = TMath::Abs(10000.*avDrft);
+ driftPath = sddLength-driftPath;
+ if(driftPath < 0) {
+ AliDebug(1, // this should be fixed at geometry level
+ Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
+ driftPath,sddLength,avDrft,dxL[0],xL[0]));
+ continue;
+ } // end if driftPath < 0
+ drTime = driftPath/driftSpeed; // drift time for segment.
+ // Sigma along the anodes for track segment.
+ sigA = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
+ sigT = sigA/driftSpeed;
+
+ drTime+=tof; // take into account Time Of Flight from production point
+ drTime-=trigDelay;
+ drTime+=timeZero;
+ timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1.001); // time bin in range 1-256 !!!
+ if(zAnode>nofAnodes) zAnode-=nofAnodes; // to have the anode number between 0. and 256.
+ iAnode = (Int_t) (1.001+zAnode); // iAnode in range 1-256 !!!!
+
+ // Peak amplitude in nanoAmpere
+ amplitude = fScaleSize*160.*depEnergy/
+ (timeStep*eVpairs*2.*acos(-1.));
+ chargeloss = 1.-cHloss*driftPath/1000.;
+ amplitude *= chargeloss;
+ amplitude *= adcscale;
+ width = 2.*nsigma/(nlookups-1);
+ // Spread the charge
+ nsplitAn = 4;
+ nsplitTb=4;
+ aStep = anodePitch/(nsplitAn*sigA);
+ aConst = zAnode*anodePitch/sigA;
+ tStep = timeStep/(nsplitTb*fScaleSize*sigT);
+ tConst = drTime/sigT;
+ // Define SDD window corresponding to the hit
+ anodeWindow = (Int_t)(nsigma*sigA/anodePitch+1);
+ timeWindow = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
+ jamin = (iAnode - anodeWindow - 2)*nsplitAn+1;
+ if(jamin <= 0) jamin = 1;
+ if(jamin > nofAnodes*nsplitAn){
+ AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode min=%d",jamin));
+ continue;