X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSsimulationSDD.cxx;h=9767f58c14bfe154804fab4d245ff135fdfd796a;hb=c77b3dc7e4eb7e6c571828c984874471eb4e3157;hp=d632d19faad637da9f4aaaf0600016526f27ca2c;hpb=13a2b50da4d399b6c7404081ea55205de2fa2bc7;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSsimulationSDD.cxx b/ITS/AliITSsimulationSDD.cxx index d632d19faad..9767f58c14b 100644 --- a/ITS/AliITSsimulationSDD.cxx +++ b/ITS/AliITSsimulationSDD.cxx @@ -18,7 +18,7 @@ #include #include #include -#include +#include #include #include @@ -32,10 +32,10 @@ #include "AliITSdigitSPD.h" #include "AliITSetfSDD.h" #include "AliITSmodule.h" +#include "AliITShit.h" #include "AliITSpList.h" -#include "AliITSresponseSDD.h" #include "AliITSCalibrationSDD.h" -#include "AliITSsegmentationSDD.h" +#include "AliITSresponseSDD.h" #include "AliITSsimulationSDD.h" #include "AliLog.h" #include "AliRun.h" @@ -49,72 +49,6 @@ ClassImp(AliITSsimulationSDD) // AliITSsimulationSDD is the simulation of SDDs. // //////////////////////////////////////////////////////////////////////// -//______________________________________________________________________ -Int_t power(Int_t b, Int_t e) { - // compute b to the e power, where both b and e are Int_ts. - Int_t power = 1,i; - - for(i=0; iGetSamples(); - Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5); - Int_t m1 = samples; - Int_t m = samples/2; - Int_t m2 = samples/m1; - Int_t i,j,k; - for(i=1; i<=l; i++) { - for(j=0; jGetWeightReal(p); - Double_t wsi = alisddetf->GetWeightImag(p); - if(direction == -1) wsi = -wsi; - Double_t xr = *(real+k+m); - Double_t xi = *(imag+k+m); - *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi); - *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr); - *(real+k) += xr; - *(imag+k) += xi; - p += m2; - } // end for k - } // end for j - m1 = m; - m /= 2; - m2 += m2; - } // end for i - - for(j=0; j= j) { - Double_t xr = *(real+j); - Double_t xi = *(imag+j); - *(real+j) = *(real+p); - *(imag+j) = *(imag+p); - *(real+p) = xr; - *(imag+p) = xi; - } // end if p>=j - } // end for j - if(direction == -1) { - for(i=0; iGetResponse(1); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); fpList = new AliITSpList( seg->Npz(), fScaleSize*seg->Npx() ); fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1); @@ -263,32 +166,17 @@ void AliITSsimulationSDD::Init(){ Double_t timeStep = (Double_t)seg->Dpx(0); if(anodePitch*(fNofMaps/2) > sddWidth) { - Warning("AliITSsimulationSDD", - "Too many anodes %d or too big pitch %f \n", - fNofMaps/2,anodePitch); + AliWarning(Form("Too many anodes %d or too big pitch %f ", + fNofMaps/2,anodePitch)); } // end if fElectronics = new AliITSetfSDD(timeStep/fScaleSize, - res->Electronics()); + simpar->GetSDDElectronics()); - char opt1[20], opt2[20]; - res->ParamOptions(opt1,opt2); - fParam = opt2; - const char *kopt=res->ZeroSuppOption(); - fD.Set(fNofMaps); - fT1.Set(fNofMaps); - fT2.Set(fNofMaps); - fTol.Set(fNofMaps); - - Bool_t write = res->OutputOption(); - if(write && strstr(kopt,"2D")) MakeTreeB(); - fITS = (AliITS*)gAlice->GetModule("ITS"); - Int_t size = fNofMaps*fMaxNofSamples; - fStream = new AliITSInStream(size); - + fInZR = new Double_t [fScaleSize*fMaxNofSamples]; fInZI = new Double_t [fScaleSize*fMaxNofSamples]; fOutZR = new Double_t [fScaleSize*fMaxNofSamples]; @@ -301,7 +189,6 @@ AliITSsimulationSDD::~AliITSsimulationSDD() { // delete fpList; delete fHitSigMap2; delete fHitNoiMap2; - delete fStream; delete fElectronics; fITS = 0; @@ -310,7 +197,6 @@ AliITSsimulationSDD::~AliITSsimulationSDD() { fHis->Delete(); delete fHis; } // end if fHis - if(fTreeB) delete fTreeB; if(fInZR) delete [] fInZR; if(fInZI) delete [] fInZI; if(fOutZR) delete [] fOutZR; @@ -332,6 +218,64 @@ void AliITSsimulationSDD::ClearMaps() { fHitSigMap2->ClearMap(); fHitNoiMap2->ClearMap(); } +//______________________________________________________________________ +void AliITSsimulationSDD::FastFourierTransform(Double_t *real, + Double_t *imag,Int_t direction) { + // Do a Fast Fourier Transform + + Int_t samples = fElectronics->GetSamples(); + Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5); + Int_t m1 = samples; + Int_t m = samples/2; + Int_t m2 = samples/m1; + Int_t i,j,k; + for(i=1; i<=l; i++) { + for(j=0; jGetWeightReal(p); + Double_t wsi = fElectronics->GetWeightImag(p); + if(direction == -1) wsi = -wsi; + Double_t xr = *(real+k+m); + Double_t xi = *(imag+k+m); + *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi); + *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr); + *(real+k) += xr; + *(imag+k) += xi; + p += m2; + } // end for k + } // end for j + m1 = m; + m /= 2; + m2 += m2; + } // end for i + for(j=0; j= j) { + Double_t xr = *(real+j); + Double_t xi = *(imag+j); + *(real+j) = *(real+p); + *(imag+j) = *(imag+p); + *(real+p) = xr; + *(imag+p) = xi; + } // end if p>=j + } // end for j + if(direction == -1) { + for(i=0; iGetResponse(1); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); Int_t nItems = pItemArray->GetEntries(); - Double_t maxadc = res->MaxAdc(); + Double_t maxadc = simpar->GetSDDMaxAdc(); Bool_t sig = kFALSE; // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl; @@ -398,14 +342,7 @@ void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){ Int_t nhits = fHits->GetEntriesFast(); InitSimulationModule( md, ev ); - - if( !nhits && fCheckNoise ) { - ChargeToSignal( fModule,kTRUE,kFALSE ); // process noise - GetNoise(); - ClearMaps(); - return; - } else - if( !nhits ) return; + if( !nhits ) return; HitsToAnalogDigits( mod ); ChargeToSignal( fModule,kTRUE,kTRUE ); // process signal + noise @@ -435,320 +372,309 @@ void AliITSsimulationSDD::FinishDigits() { if( fCrosstalkFlag ) ApplyCrosstalk(fModule); AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - const char *kopt = res->GetZeroSuppOption(); - ZeroSuppression( kopt ); + Bool_t isZeroSupp = res->GetZeroSupp(); + if (isZeroSupp) Compress2D(); + else StoreAllDigits(); } //______________________________________________________________________ void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) { // create maps to build the lists of tracks for each digit AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1); AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); TObjArray *hits = mod->GetHits(); - Int_t nhits = hits->GetEntriesFast(); - - // Int_t arg[6] = {0,0,0,0,0,0}; - Int_t nofAnodes = fNofMaps/2; - Double_t sddLength = seg->Dx(); - Double_t sddWidth = seg->Dz(); - Double_t anodePitch = seg->Dpz(0); - Double_t timeStep = seg->Dpx(0); - Double_t driftSpeed ; // drift velocity (anode dependent) - //Float_t maxadc = res->GetMaxAdc(); - //Float_t topValue = res->GetDynamicRange(); - Double_t norm = res->GetMaxAdc()/res->GetDynamicRange(); // maxadc/topValue; - Double_t cHloss = res->GetChargeLoss(); - Float_t dfCoeff, s1; res->DiffCoeff(dfCoeff,s1); // Signal 2d Shape - Double_t eVpairs = res->GetGeVToCharge()*1.0E9; // 3.6 eV by def. - Double_t nsigma = res->GetNSigmaIntegration(); // - Int_t nlookups = res->GetGausNLookUp(); // - Float_t jitter = res->GetJitterError(); // - - // Piergiorgio's part (apart for few variables which I made float - // when i thought that can be done - // Fill detector maps with GEANT hits - // loop over hits in the module - - const Float_t kconv = 1.0e+6; // GeV->KeV - Int_t itrack = 0; - Int_t iWing; // which detector wing/side. - Int_t ii,kk,ka,kt; // loop indexs - Int_t ia,it,index; // sub-pixel integration indexies - Int_t iAnode; // anode number. - Int_t timeSample; // time buckett. - Int_t anodeWindow; // anode direction charge integration width - Int_t timeWindow; // time direction charge integration width - Int_t jamin,jamax; // anode charge integration window - Int_t jtmin,jtmax; // time charge integration window - Int_t ndiv; // Anode window division factor. - Int_t nsplit; // the number of splits in anode and time windows==1. - Int_t nOfSplits; // number of times track length is split into - Float_t nOfSplitsF; // Floating point version of nOfSplits. - Float_t kkF; // Floating point version of loop index kk. - Double_t pathInSDD; // Track length in SDD. - Double_t drPath; // average position of track in detector. in microns - Double_t drTime; // Drift time - Double_t nmul; // drift time window multiplication factor. - Double_t avDrft; // x position of path length segment in cm. - Double_t avAnode; // Anode for path length segment in Anode number (float) - Double_t zAnode; // Floating point anode number. - Double_t driftPath; // avDrft in microns. - Double_t width; // width of signal at anodes. - Double_t depEnergy; // Energy deposited in this GEANT step. - Double_t xL[3],dxL[3]; // local hit coordinates and diff. - Double_t sigA; // sigma of signal at anode. - Double_t sigT; // sigma in time/drift direction for track segment - Double_t aStep,aConst; // sub-pixel size and offset anode - Double_t tStep,tConst; // sub-pixel size and offset time - Double_t amplitude; // signal amplitude for track segment in nanoAmpere - Double_t chargeloss; // charge loss for track segment. - Double_t anodeAmplitude; // signal amplitude in anode direction - Double_t aExpo; // exponent of Gaussian anode direction - Double_t timeAmplitude; // signal amplitude in time direction - Double_t tExpo; // exponent of Gaussian time direction - // Double_t tof; // Time of flight in ns of this step. - - for(ii=0; iiLineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2], - depEnergy,itrack)) continue; - Float_t xloc=xL[0]; - if(xloc>0) iWing=0; // left side, carlos channel 0 - else iWing=1; // right side - - Float_t zloc=xL[2]+0.5*dxL[2]; - zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511. - driftSpeed = res->GetDriftSpeedAtAnode(zAnode); - if(timeStep*fMaxNofSamples < sddLength/driftSpeed) { - AliWarning("Time Interval > Allowed Time Interval\n"); - } - depEnergy *= kconv; - - // scale path to simulate a perpendicular track + Int_t nhits = hits->GetEntriesFast(); + + // Int_t arg[6] = {0,0,0,0,0,0}; + Int_t nofAnodes = fNofMaps/2; + Double_t sddLength = seg->Dx(); + Double_t anodePitch = seg->Dpz(0); + Double_t timeStep = seg->Dpx(0); + Double_t driftSpeed ; // drift velocity (anode dependent) + Double_t nanoampToADC = simpar->GetSDDMaxAdc()/simpar->GetSDDDynamicRange(); // maxadc/topValue; + Double_t cHloss = simpar->GetSDDChargeLoss(); + Float_t dfCoeff, s1; + simpar->GetSDDDiffCoeff(dfCoeff,s1); // Signal 2d Shape + Double_t eVpairs = simpar->GetGeVToCharge()*1.0E9; // 3.6 eV by def. + Double_t nsigma = simpar->GetNSigmaIntegration(); // + Int_t nlookups = simpar->GetGausNLookUp(); // + Float_t jitter = simpar->GetSDDJitterError(); // + Float_t trigDelay = simpar->GetSDDTrigDelay(); // compensation for MC time zero + Float_t timeZero=fDetType->GetResponseSDD()->GetTimeZero(fModule); + + // Piergiorgio's part (apart for few variables which I made float + // when i thought that can be done + // Fill detector maps with GEANT hits + // loop over hits in the module + + const Float_t kconv = 1.0e+6; // GeV->KeV + Int_t itrack = 0; + Int_t iWing; // which detector wing/side. + Int_t ii,kk,ka,kt; // loop indexs + Int_t ia,it,index; // sub-pixel integration indexies + Int_t iAnode; // anode number. + Int_t timeSample; // time buckett. + Int_t anodeWindow; // anode direction charge integration width + Int_t timeWindow; // time direction charge integration width + Int_t jamin,jamax; // anode charge integration window + Int_t jtmin,jtmax; // time charge integration window + Int_t nsplitAn; // the number of splits in anode and time windows + Int_t nsplitTb; // the number of splits in anode and time windows + Int_t nOfSplits; // number of times track length is split into + Float_t nOfSplitsF; // Floating point version of nOfSplits. + Float_t kkF; // Floating point version of loop index kk. + Double_t pathInSDD; // Track length in SDD. + Double_t drPath; // average position of track in detector. in microns + Double_t drTime; // Drift time + Double_t avDrft; // x position of path length segment in cm. + Double_t avAnode; // Anode for path length segment in Anode number (float) + Double_t zAnode; // Floating point anode number. + Double_t driftPath; // avDrft in microns. + Double_t width; // width of signal at anodes. + Double_t depEnergy; // Energy deposited in this GEANT step. + Double_t xL[3],dxL[3]; // local hit coordinates and diff. + Double_t sigA; // sigma of signal at anode. + Double_t sigT; // sigma in time/drift direction for track segment + Double_t aStep,aConst; // sub-pixel size and offset anode + Double_t tStep,tConst; // sub-pixel size and offset time + Double_t amplitude; // signal amplitude for track segment in nanoAmpere + Double_t chargeloss; // charge loss for track segment. + Double_t anodeAmplitude; // signal amplitude in anode direction + Double_t aExpo; // exponent of Gaussian anode direction + Double_t timeAmplitude; // signal amplitude in time direction + Double_t tExpo; // exponent of Gaussian time direction + Double_t tof; // Time of flight in ns of this step. + + for(ii=0; iiLineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2], + depEnergy,itrack)) continue; + Float_t xloc=xL[0]; + if(xloc>0) iWing=0; // left side, carlos channel 0 + else iWing=1; // right side + + Float_t zloc=xL[2]+0.5*dxL[2]; + zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511. + driftSpeed = res->GetDriftSpeedAtAnode(zAnode); + if(timeStep*fMaxNofSamples < sddLength/driftSpeed) { + AliWarning("Time Interval > Allowed Time Interval"); + } + depEnergy *= kconv; + if (!depEnergy) { + AliDebug(1, + Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!", + itrack,ii,mod->GetIndex())); + continue; // continue if the particle did not lose energy // passing through detector - if (!depEnergy) { - AliDebug(1, - Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!", - itrack,ii,mod->GetIndex())); - continue; - } // end if !depEnergy - - xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); // - pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]); + } // end if !depEnergy + + tof=0.; + AliITShit* h=(AliITShit*)hits->At(ii); + if(h){ + tof=h->GetTOF()*1E9; + AliDebug(1,Form("TOF for hit %d on mod %d (particle %d)=%g",ii,fModule,h->Track(),tof)); + } - if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); } - drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5); - drPath = sddLength-drPath; - if(drPath < 0) { + xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); // + pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]); + + if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); } + drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5); + drPath = sddLength-drPath; + if(drPath < 0) { + AliDebug(1, // this should be fixed at geometry level + Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e", + drPath,sddLength,dxL[0],xL[0])); + continue; + } // end if drPath < 0 + + // Compute number of segments to brake step path into + drTime = drPath/driftSpeed; // Drift Time + sigA = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes + // calcuate the number of time the path length should be split into. + nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA); + if(fFlag) nOfSplits = 1; + + // loop over path segments, init. some variables. + depEnergy /= nOfSplits; + nOfSplitsF = (Float_t) nOfSplits; + Float_t theAverage=0.,theSteps=0.; + for(kk=0;kkGetAnodeFromLocal(avDrft,avAnode); + driftSpeed = res->GetDriftSpeedAtAnode(zAnode); + driftPath = TMath::Abs(10000.*avDrft); + driftPath = sddLength-driftPath; + if(driftPath < 0) { AliDebug(1, // this should be fixed at geometry level - Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e", - drPath,sddLength,dxL[0],xL[0])); + 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 drPath < 0 - - // Compute number of segments to brake step path into - drTime = drPath/driftSpeed; // Drift Time - sigA = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes - // calcuate the number of time the path length should be split into. - nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA); - if(fFlag) nOfSplits = 1; + } // 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 !!!! - // loop over path segments, init. some variables. - depEnergy /= nOfSplits; - nOfSplitsF = (Float_t) nOfSplits; - Float_t theAverage=0.,theSteps=0.; - for(kk=0;kkGetAnodeFromLocal(avDrft,avAnode); - driftSpeed = res->GetDriftSpeedAtAnode(zAnode); - 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. - timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1); // time bin in range 1-256 !!! - if(timeSample > fScaleSize*fMaxNofSamples) { - AliWarning(Form("Wrong Time Sample: %e",timeSample)); - continue; - } // end if timeSample > fScaleSize*fMaxNoofSamples - - if(zAnode>nofAnodes) zAnode-=nofAnodes; // to have the anode number between 0. and 256. - if(zAnode*anodePitch > sddWidth || zAnode*anodePitch < 0.) - AliWarning(Form("Exceeding sddWidth=%e Z = %e",sddWidth,zAnode*anodePitch)); - iAnode = (Int_t) (1.+zAnode); // iAnode in range 1-256 !!!! - if(iAnode < 1 || iAnode > nofAnodes) { - AliWarning(Form("Wrong iAnode: 1<%d>%d (xanode=%e)",iAnode,nofAnodes, zAnode)); - continue; - } // end if iAnode < 1 || iAnode > nofAnodes - - // 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) : - - // Sigma along the anodes for track segment. - sigA = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1); - sigT = sigA/driftSpeed; // Peak amplitude in nanoAmpere - amplitude = fScaleSize*160.*depEnergy/ - (timeStep*eVpairs*2.*acos(-1.)*sigT*sigA); - amplitude *= timeStep/25.; // WARNING!!!!! Amplitude scaling to - // account for clock variations - // (reference value: 40 MHz) - chargeloss = 1.-cHloss*driftPath/1000.; - amplitude *= chargeloss; - width = 2.*nsigma/(nlookups-1); - // Spread the charge - // Pixel index - ndiv = 2; - nmul = 3.; - if(drTime > 1200.) { - ndiv = 4; - nmul = 1.5; - } // end if drTime > 1200. - // Sub-pixel index - nsplit = 4; // hard-wired //nsplit=4;nsplit = (nsplit+1)/2*2; - // Sub-pixel size see computation of aExpo and tExpo. - aStep = anodePitch/(nsplit*fScaleSize*sigA); - aConst = zAnode*anodePitch/sigA; - tStep = timeStep/(nsplit*fScaleSize*sigT); - tConst = drTime/sigT; - // Define SDD window corresponding to the hit - anodeWindow = (Int_t)(fScaleSize*nsigma*sigA/anodePitch+1); - timeWindow = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.); - jamin = (iAnode - anodeWindow/ndiv - 2)*fScaleSize*nsplit +1; - jamax = (iAnode + anodeWindow/ndiv + 1)*fScaleSize*nsplit; - if(jamin <= 0) jamin = 1; - if(jamax > fScaleSize*nofAnodes*nsplit) - jamax = fScaleSize*nofAnodes*nsplit; - // jtmin and jtmax are Hard-wired - jtmin = (Int_t)(timeSample-timeWindow*nmul-1)*nsplit+1; - jtmax = (Int_t)(timeSample+timeWindow*nmul)*nsplit; - if(jtmin <= 0) jtmin = 1; - if(jtmax > fScaleSize*fMaxNofSamples*nsplit) - jtmax = fScaleSize*fMaxNofSamples*nsplit; - // Spread the charge in the anode-time window - for(ka=jamin; ka <=jamax; ka++) { - ia = (ka-1)/(fScaleSize*nsplit) + 1; - if(ia <= 0) { - Warning("HitsToAnalogDigits","ia < 1: "); - continue; - } // end if - if(ia > nofAnodes) ia = nofAnodes; - aExpo = (aStep*(ka-0.5)-aConst); - if(TMath::Abs(aExpo) > nsigma) anodeAmplitude = 0.; - else { - Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5); - anodeAmplitude = amplitude*res->GetGausLookUp(theBin); - } // end if TMath::Abs(aEspo) > nsigma - // index starts from 0 - index = iWing*nofAnodes+ia-1; - if(anodeAmplitude){ - for(kt=jtmin; kt<=jtmax; kt++) { - it = (kt-1)/nsplit+1; // it starts from 1 - if(it<=0){ - Warning("HitsToAnalogDigits","it < 1:"); - continue; - } // end if - if(it>fScaleSize*fMaxNofSamples) - it = fScaleSize*fMaxNofSamples; - tExpo = (tStep*(kt-0.5)-tConst); - if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.; - else { - Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5); - timeAmplitude = anodeAmplitude*res->GetGausLookUp(theBin); - } // end if TMath::Abs(tExpo) > nsigma - // build the list of Sdigits for this module - // arg[0] = index; - // arg[1] = it; - // arg[2] = itrack; // track number - // arg[3] = ii-1; // hit number. - timeAmplitude *= norm; - timeAmplitude *= 10; - // ListOfFiredCells(arg,timeAmplitude,alst,padr); - Double_t charge = timeAmplitude; - charge += fHitMap2->GetSignal(index,it-1); - fHitMap2->SetHit(index, it-1, charge); - fpList->AddSignal(index,it-1,itrack,ii-1, - mod->GetIndex(),timeAmplitude); - fAnodeFire[index] = kTRUE; - } // end loop over time in window - } // end if anodeAmplitude - } // loop over anodes in window - } // end loop over "sub-hits" - } // end loop over hits + amplitude = fScaleSize*160.*depEnergy/ + (timeStep*eVpairs*2.*acos(-1.)); + chargeloss = 1.-cHloss*driftPath/1000.; + amplitude *= chargeloss; + 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; + } + jamax = (iAnode + anodeWindow + 2)*nsplitAn; + if(jamax > nofAnodes*nsplitAn) jamax = nofAnodes*nsplitAn; + if(jamax <=0){ + AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode max=%d",jamax)); + continue; + } + jtmin = (Int_t)(timeSample-timeWindow-2)*nsplitTb+1; + if(jtmin <= 0) jtmin = 1; + if(jtmin > fScaleSize*fMaxNofSamples*nsplitTb){ + AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample min=%d tof=%f",jtmin,tof)); + continue; + } + jtmax = (Int_t)(timeSample+timeWindow+2)*nsplitTb; + if(jtmax > fScaleSize*fMaxNofSamples*nsplitTb) jtmax = fScaleSize*fMaxNofSamples*nsplitTb; + if(jtmax <= 0){ + AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample max=%d tof=%f",jtmax,tof)); + continue; + } + + // Spread the charge in the anode-time window + for(ka=jamin; ka <=jamax; ka++) { + ia = (ka-1)/nsplitAn + 1; + if(ia <= 0) ia=1; + if(ia > nofAnodes) ia = nofAnodes; + aExpo = (aStep*(ka-0.5)-aConst); + if(TMath::Abs(aExpo) > nsigma) anodeAmplitude = 0.; + else { + Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5); + anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin); + } + // index starts from 0 + index = iWing*nofAnodes+ia-1; + if(anodeAmplitude){ + for(kt=jtmin; kt<=jtmax; kt++) { + it = (kt-1)/nsplitTb+1; // it starts from 1 + if(it<=0) it=1; + if(it>fScaleSize*fMaxNofSamples) + it = fScaleSize*fMaxNofSamples; + tExpo = (tStep*(kt-0.5)-tConst); + if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.; + else { + Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5); + timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin)*aStep*tStep; + } + timeAmplitude *= nanoampToADC; + // ListOfFiredCells(arg,timeAmplitude,alst,padr); + Double_t charge = timeAmplitude; + charge += fHitMap2->GetSignal(index,it-1); + fHitMap2->SetHit(index, it-1, charge); + fpList->AddSignal(index,it-1,itrack,ii-1, + mod->GetIndex(),timeAmplitude); + fAnodeFire[index] = kTRUE; + } // end loop over time in window + } // end if anodeAmplitude + } // loop over anodes in window + } // end loop over "sub-hits" + } // end loop over hits } //____________________________________________ -void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signal ) { - // Adds a Digit. - Int_t size = AliITSdigit::GetNTracks(); - - Int_t digits[3]; - Int_t * tracks = new Int_t[size]; - Int_t * hits = new Int_t[size]; - Float_t phys; - Float_t * charges = new Float_t[size]; - - digits[0] = i; - digits[1] = j; - digits[2] = signal; - - AliITSpListItem *pItem = fpList->GetpListItem( i, j ); - if( pItem == 0 ) { - phys = 0.0; - for( Int_t l=0; lGetTrack( 0 ); - if( idtrack >= 0 ) phys = pItem->GetSignal(); - else phys = 0.0; - - for( Int_t l=0; lGetMaxKept()) { - tracks[l] = pItem->GetTrack( l ); - hits[l] = pItem->GetHit( l ); - charges[l] = pItem->GetSignal( l ); - }else{ - tracks[l] = -3; - hits[l] = -1; - charges[l] = 0.0; - }// end for if +void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signalc, Int_t signale) { + // Adds a Digit. + Int_t size = AliITSdigit::GetNTracks(); + + Int_t digits[3]; + Int_t * tracks = new Int_t[size]; + Int_t * hits = new Int_t[size]; + Float_t phys; + Float_t * charges = new Float_t[size]; + + digits[0] = i; + digits[1] = j; + digits[2] = signalc; + + AliITSpListItem *pItem = fpList->GetpListItem( i, j ); + if( pItem == 0 ) { + phys = 0.0; + for( Int_t l=0; lGetTrack( 0 ); + if( idtrack >= 0 ) phys = pItem->GetSignal(); + else phys = 0.0; + + for( Int_t l=0; lGetMaxKept()) { + tracks[l] = pItem->GetTrack( l ); + hits[l] = pItem->GetHit( l ); + charges[l] = pItem->GetSignal( l ); + }else{ + tracks[l] = -3; + hits[l] = -1; + charges[l] = 0.0; + }// end for if + } - fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges ); - delete [] tracks; - delete [] hits; - delete [] charges; + fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges, signale ); + delete [] tracks; + delete [] hits; + delete [] charges; } //______________________________________________________________________ void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAddGain) { // add baseline, noise, gain, electronics and ADC saturation effects // apply dead channels - char opt1[20], opt2[20]; AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod); - res->GetParamOptions(opt1,opt2); Double_t baseline=0; Double_t noise=0; Double_t gain=0; Float_t contrib=0; Int_t i,k,kk; - Float_t maxadc = res->GetMaxAdc(); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); + Float_t maxadc = simpar->GetSDDMaxAdc(); + Int_t nGroup=fScaleSize; + if(res->IsAMAt20MHz()){ + nGroup=fScaleSize/2; + } for (i=0;iGetBaseline(i); noise = res->GetNoise(i); - gain = res->GetChannelGain(i); + gain = res->GetChannelGain(i)/fDetType->GetAverageGainSDD(); if(res->IsBad()) gain=0.; if( res->IsChipBad(res->GetChip(i)) )gain=0.; for(k=0; k= maxadc) newcont = maxadc -1; if(newcont >= baseline){ - Warning("","newcont=%d>=baseline=%d",newcont,baseline); + Warning("","newcont=%f>=baseline=%f",newcont,baseline); } // end if // back to analog: ? fHitMap2->SetHit(i,k,newcont); } // end for k }else{ - FastFourierTransform(fElectronics,&fInZR[0],&fInZI[0],1); + FastFourierTransform(&fInZR[0],&fInZI[0],1); for(k=0; kGetTraFunReal(k); Double_t iw = fElectronics->GetTraFunImag(k); fOutZR[k] = fInZR[k]*rw - fInZI[k]*iw; fOutZI[k] = fInZR[k]*iw + fInZI[k]*rw; } // end for k - FastFourierTransform(fElectronics,&fOutZR[0],&fOutZI[0],-1); + FastFourierTransform(&fOutZR[0],&fOutZI[0],-1); for(k=0; k maxcont1) maxcont1 = newcont1; } // end for kk @@ -805,17 +731,9 @@ void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAdd void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) { // function add the crosstalk effect to signal // temporal function, should be checked...!!! - AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1); - Int_t fNofMaps = seg->Npz(); - Int_t fMaxNofSamples = seg->Npx(); - // create and inizialice crosstalk map Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1]; - if( ctk == NULL ) { - Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" ); - return; - } memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) ); AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod); for( Int_t z=0; z 2 && i < fMaxNofSamples-2 ) dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 ) @@ -888,33 +801,6 @@ void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) { delete [] ctk; } -//______________________________________________________________________ -void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl, - Int_t &th) const{ - // Returns the compression alogirthm parameters - db=fD[i]; - tl=fT1[i]; - th=fT2[i]; -} -//______________________________________________________________________ -void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl) const{ - // returns the compression alogirthm parameters - - db=fD[i]; - tl=fT1[i]; - -} -//______________________________________________________________________ -void AliITSsimulationSDD::SetCompressParam(){ - // Sets the compression alogirthm parameters - AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - for(Int_t ian = 0; ianGetBaseline(ian)); - fT1[ian] = (Int_t)(2.*calibr->GetNoiseAfterElectronics(ian)+0.5); - fT2[ian] = 0; // used by 2D clustering - not defined yet - fTol[ian] = 0; // used by 2D clustering - not defined yet - } -} //______________________________________________________________________ Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const { // To the 10 to 8 bit lossive compression. @@ -926,353 +812,100 @@ Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const { if (signal < 1024) return (224+((signal-512)>>4)); return 0; } -/* -//______________________________________________________________________ -AliITSMap* AliITSsimulationSDD::HitMap(Int_t i){ - //Return the correct map. - - return ((i==0)? fHitMap1 : fHitMap2); -} -*/ //______________________________________________________________________ -void AliITSsimulationSDD::ZeroSuppression(const char *option) { - // perform the zero suppresion - if (strstr(option,"2D")) { - //Init2D(); // activate if param change module by module - Compress2D(); - } else if (strstr(option,"1D")) { - //Init1D(); // activate if param change module by module - Compress1D(); - } else StoreAllDigits(); -} -//______________________________________________________________________ -void AliITSsimulationSDD::Init2D(){ - // read in and prepare arrays: fD, fT1, fT2 - // savemu[nanodes], savesigma[nanodes] - // read baseline and noise from file - either a .root file and in this - // case data should be organised in a tree with one entry for each - // module => reading should be done accordingly - // or a classic file and do smth. like this ( code from Davide C. and - // Albert W.) : - // Read 2D zero-suppression parameters for SDD - - if (!strstr(fParam.Data(),"file")) return; +Int_t AliITSsimulationSDD::Convert8to10(Int_t signal) const { + // Decompression from 8 to 10 bit - Int_t na,pos,tempTh; - Float_t mu,sigma; - Float_t *savemu = new Float_t [fNofMaps]; - Float_t *savesigma = new Float_t [fNofMaps]; - char input[100],basel[100],par[100]; - char *filtmp; - Double_t tmp1,tmp2; - AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - - res->Thresholds(tmp1,tmp2); - Int_t minval = static_cast(tmp1); - - res->Filenames(input,basel,par); - fFileName = par; - // - filtmp = gSystem->ExpandPathName(fFileName.Data()); - FILE *param = fopen(filtmp,"r"); - na = 0; - - if(param) { - while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) { - if (pos != na+1) { - Error("Init2D","Anode number not in increasing order!",filtmp); - exit(1); - } // end if pos != na+1 - savemu[na] = mu; - savesigma[na] = sigma; - if ((2.*sigma) < mu) { - fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5); - mu = 2.0 * sigma; - } else fD[na] = 0; - tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval; - if (tempTh < 0) tempTh=0; - fT1[na] = tempTh; - tempTh = (Int_t)floor(mu+3.0*sigma+0.5) - minval; - if (tempTh < 0) tempTh=0; - fT2[na] = tempTh; - na++; - } // end while - } else { - Error("Init2D","THE FILE %s DOES NOT EXIST !",filtmp); - exit(1); - } // end if(param) + if (signal < 0 || signal > 255) { + AliWarning(Form("Signal value %d out of range",signal)); + return 0; + } // end if signal <0 || signal >255 - fclose(param); - delete [] filtmp; - delete [] savemu; - delete [] savesigma; + if (signal < 128) return signal; + if (signal < 192) { + if (TMath::Odd(signal)) return (128+((signal-128)<<1)); + else return (128+((signal-128)<<1)+1); + } // end if signal < 192 + if (signal < 224) { + if (TMath::Odd(signal)) return (256+((signal-192)<<3)+3); + else return (256+((signal-192)<<3)+4); + } // end if signal < 224 + if (TMath::Odd(signal)) return (512+((signal-224)<<4)+7); + return (512+((signal-224)<<4)+8); } //______________________________________________________________________ void AliITSsimulationSDD::Compress2D(){ - // simple ITS cluster finder -- online zero-suppression conditions - - Int_t db,tl,th; - Double_t tmp1,tmp2; - AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - - res->Thresholds(tmp1,tmp2); - Int_t minval = static_cast(tmp1); - Bool_t write = res->OutputOption(); - Bool_t do10to8 = res->Do10to8(); - Int_t nz, nl, nh, low, i, j; - SetCompressParam(); - for (i=0; iGetSignal(i,j)); - signal -= db; // if baseline eq. is done here - if (signal <= 0) {nz++; continue;} - if ((signal - tl) < minval) low++; - if ((signal - th) >= minval) { - nh++; - Bool_t cond=kTRUE; - FindCluster(i,j,signal,minval,cond); - if(cond && j && - ((TMath::Abs(fHitMap2->GetSignal(i,j-1))-th)>=minval)){ - if(do10to8) signal = Convert10to8(signal); - AddDigit(i,j,signal); - } // end if cond&&j&&() - } else if ((signal - tl) >= minval) nl++; - } // end for j loop time samples - if (write) TreeB()->Fill(nz,nl,nh,low,i+1); - } //end for i loop anodes - - char hname[30]; - if (write) { - sprintf(hname,"TNtuple%d_%d",fModule,fEvent); - TreeB()->Write(hname); - // reset tree - TreeB()->Reset(); - } // end if write -} -//______________________________________________________________________ -void AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal, - Int_t minval,Bool_t &cond){ - // Find clusters according to the online 2D zero-suppression algorithm - AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1); - - Bool_t do10to8 = res->Do10to8(); - Bool_t high = kFALSE; - - fHitMap2->FlagHit(i,j); - // - // check the online zero-suppression conditions - // - const Int_t kMaxNeighbours = 4; - Int_t nn; - Int_t dbx,tlx,thx; - Int_t xList[kMaxNeighbours], yList[kMaxNeighbours]; - seg->Neighbours(i,j,&nn,xList,yList); - Int_t in,ix,iy,qns; - for (in=0; inTestHit(ix,iy)==kUnused) { - CompressionParam(ix,dbx,tlx,thx); - Int_t qn = (Int_t)(fHitMap2->GetSignal(ix,iy)); - qn -= dbx; // if baseline eq. is done here - if ((qn-tlx) < minval) { - fHitMap2->FlagHit(ix,iy); - continue; - } else { - if ((qn - thx) >= minval) high=kTRUE; - if (cond) { - if(do10to8) signal = Convert10to8(signal); - AddDigit(i,j,signal); - } // end if cond - if(do10to8) qns = Convert10to8(qn); - else qns=qn; - if (!high) AddDigit(ix,iy,qns); - cond=kFALSE; - if(!high) fHitMap2->FlagHit(ix,iy); - } // end if qn-tlx < minval - } // end if TestHit - } // end for in loop over neighbours + // 2D zero-suppression algorithm as described in ALICE-INT-1999-28 V10 + AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); + for (Int_t iWing=0; iWing<2; iWing++) { + Int_t tL=res->GetZSLowThreshold(iWing); + Int_t tH=res->GetZSHighThreshold(iWing); + for (Int_t i=0; iGetSignal(ian,itb); + if(cC<=tL) continue; + nLow++; // cC is greater than tL + if(cC>tH) nHigh++; + // N + // Get "quintuple": WCE + // S + Float_t wW=0.; + if(itb>0) wW=fHitMap2->GetSignal(ian,itb-1); + if(wW>tL) nLow++; + if(wW>tH) nHigh++; + Float_t eE=0.; + if(itbGetSignal(ian,itb+1); + if(eE>tL) nLow++; + if(eE>tH) nHigh++; + Float_t nN=0.; + if(i<(fNofMaps/2-1)) nN=fHitMap2->GetSignal(ian+1,itb); + if(nN>tL) nLow++; + if(nN>tH) nHigh++; + Float_t sS=0.; + if(i>0) sS=fHitMap2->GetSignal(ian-1,itb); + if(sS>tL) nLow++; + if(sS>tH) nHigh++; + + if(nLow>=2 && nHigh>=1){ + Int_t signal=(Int_t)cC; + Int_t signalc = Convert10to8(signal); + Int_t signale = Convert8to10(signalc); + signalc-=tL; // subtract low threshold after 10 to 8 bit compression + if(signalc>=4) AddDigit(ian,itb,signalc,signale); // store C + } + } + } + } } -//______________________________________________________________________ -void AliITSsimulationSDD::Init1D(){ - // this is just a copy-paste of input taken from 2D algo - // Torino people should give input - // Read 1D zero-suppression parameters for SDD - - if (!strstr(fParam.Data(),"file")) return; - - Int_t na,pos,tempTh; - Float_t mu,sigma; - Float_t *savemu = new Float_t [fNofMaps]; - Float_t *savesigma = new Float_t [fNofMaps]; - char input[100],basel[100],par[100]; - char *filtmp; - Double_t tmp1,tmp2; - AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - res->Thresholds(tmp1,tmp2); - Int_t minval = static_cast(tmp1); - res->Filenames(input,basel,par); - fFileName=par; - - // set first the disable and tol param - SetCompressParam(); - // - filtmp = gSystem->ExpandPathName(fFileName.Data()); - FILE *param = fopen(filtmp,"r"); - na = 0; - - if (param) { - fscanf(param,"%d %d %d %d ", &fT2[0], &fT2[1], &fTol[0], &fTol[1]); - while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) { - if (pos != na+1) { - Error("Init1D","Anode number not in increasing order!",filtmp); - exit(1); - } // end if pos != na+1 - savemu[na]=mu; - savesigma[na]=sigma; - if ((2.*sigma) < mu) { - fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5); - mu = 2.0 * sigma; - } else fD[na] = 0; - tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval; - if (tempTh < 0) tempTh=0; - fT1[na] = tempTh; - na++; - } // end while - } else { - Error("Init1D","THE FILE %s DOES NOT EXIST !",filtmp); - exit(1); - } // end if(param) - - fclose(param); - delete [] filtmp; - delete [] savemu; - delete [] savesigma; -} -//______________________________________________________________________ -void AliITSsimulationSDD::Compress1D(){ - // 1D zero-suppression algorithm (from Gianluca A.) - Int_t dis,tol,thres,decr,diff; - UChar_t *str=fStream->Stream(); - Int_t counter=0; - AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - - Bool_t do10to8=res->Do10to8(); - Int_t last=0; - Int_t k,i,j; - SetCompressParam(); - for (k=0; k<2; k++) { - tol = Tolerance(k); - dis = Disable(k); - for (i=0; iGetSignal(idx,j)); - signal -= decr; // if baseline eq. - if(do10to8) signal = Convert10to8(signal); - if (signal <= thres) { - signal=0; - diff=128; - last=0; - // write diff in the buffer for HuffT - str[counter]=(UChar_t)diff; - counter++; - continue; - } // end if signal <= thres - diff=signal-last; - if (diff > 127) diff=127; - if (diff < -128) diff=-128; - if (signal < dis) { - // tol has changed to 8 possible cases ? - one can write - // this if(TMath::Abs(diff)CheckCount(counter); - - // open file and write out the stream of diff's - static Bool_t open=kTRUE; - static TFile *outFile; - Bool_t write = res->OutputOption(); - TDirectory *savedir = gDirectory; - - if (write ) { - if(open) { - SetFileName("stream.root"); - cout<<"filename "<cd(); - fStream->Write(); - } // endif write - - fStream->ClearStream(); - - // back to galice.root file - if(savedir) savedir->cd(); -} //______________________________________________________________________ void AliITSsimulationSDD::StoreAllDigits(){ - // if non-zero-suppressed data - AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - - Bool_t do10to8 = res->Do10to8(); - Int_t i, j, digits[3]; - - for (i=0; iGetSignal(i,j)); - if(do10to8) signal = Convert10to8(signal); - digits[0] = i; - digits[1] = j; - digits[2] = signal; - fITS->AddRealDigit(1,digits); - } // end for j - } // end for i + // store digits for non-zero-suppressed data + for (Int_t ian=0; ianGetSignal(ian,itb)); + Int_t signalc = Convert10to8(signal); + Int_t signale = Convert8to10(signalc); + AddDigit(ian,itb,signalc,signale); + } + } } //______________________________________________________________________ void AliITSsimulationSDD::CreateHistograms(Int_t scale){ - // Creates histograms of maps for debugging - Int_t i; - - fHis=new TObjArray(fNofMaps); - for (i=0;iAddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples, - 0.,(Float_t) scale*fMaxNofSamples), i); - } // end for i + // Creates histograms of maps for debugging + Int_t i; + + fHis=new TObjArray(fNofMaps); + for (i=0;iAddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples, + 0.,(Float_t) scale*fMaxNofSamples), i); + } // end for i } //______________________________________________________________________ void AliITSsimulationSDD::FillHistograms(){ @@ -1328,57 +961,6 @@ void AliITSsimulationSDD::WriteToFile(TFile *hfile) { return; } //______________________________________________________________________ -Float_t AliITSsimulationSDD::GetNoise() { - // Returns the noise value - //Bool_t do10to8=GetResp()->Do10to8(); - //noise will always be in the liniar part of the signal - Int_t decr; - Int_t threshold = fT1[0]; - char opt1[20], opt2[20]; - AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule); - SetCompressParam(); - res->GetParamOptions(opt1,opt2); - fParam=opt2; - Double_t noise,baseline; - //GetBaseline(fModule); - - TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2"); - if(c2) delete c2->GetPrimitive("noisehist"); - if(c2) delete c2->GetPrimitive("anode"); - else c2=new TCanvas("c2"); - c2->cd(); - c2->SetFillColor(0); - - TH1F *noisehist = new TH1F("noisehist","noise",100,0.,(float)2*threshold); - TH1F *anode = new TH1F("anode","Anode Projection",fMaxNofSamples,0., - (float)fMaxNofSamples); - Int_t i,k; - for (i=0;iGetBaseline(i); - noise = res->GetNoise(i); - anode->Reset(); - for (k=0;kGetSignal(i,k); - //if (signal <= (float)threshold) noisehist->Fill(signal-baseline); - if (signal <= (float)(threshold+decr)) noisehist->Fill(signal); - anode->Fill((float)k,signal); - } // end for k - anode->Draw(); - c2->Update(); - } // end for i - TF1 *gnoise = new TF1("gnoise","gaus",0.,threshold); - noisehist->Fit("gnoise","RQ"); - noisehist->Draw(); - c2->Update(); - Float_t mnoise = gnoise->GetParameter(1); - cout << "mnoise : " << mnoise << endl; - Float_t rnoise = gnoise->GetParameter(2); - cout << "rnoise : " << rnoise << endl; - delete noisehist; - return rnoise; -} -//______________________________________________________________________ void AliITSsimulationSDD::WriteSDigits(){ // Fills the Summable digits Tree static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS"); @@ -1413,9 +995,8 @@ void AliITSsimulationSDD::PrintStatus() const { cout << " Silicon Drift Detector Simulation Parameters " << endl; cout << "**************************************************" << endl; cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl; - cout << "Flag for noise checking: " << (Int_t) fCheckNoise << endl; cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl; - cout << "Number pf Anodes used: " << fNofMaps << endl; + cout << "Number of Anodes used: " << fNofMaps << endl; cout << "Number of Time Samples: " << fMaxNofSamples << endl; cout << "Scale size factor: " << fScaleSize << endl; cout << "**************************************************" << endl;