/*
$Log$
+Revision 1.53 2001/05/31 18:52:24 barbera
+Bari model becomes the default
+
+Revision 1.53 2001/05/30 07:52:24 hristov
+TPC and CONTAINERS included in the search path
+
Revision 1.52 2001/05/30 06:04:58 hristov
Changes made to be consitant with changes in TPC tracking classes (B.Nilsen)
#include "AliITSresponse.h"
#include "AliITSsegmentationSPD.h"
#include "AliITSresponseSPD.h"
-#include "AliITSresponseSPDbari.h"
+#include "AliITSresponseSPDdubna.h"
#include "AliITSsegmentationSDD.h"
#include "AliITSresponseSDD.h"
#include "AliITSsegmentationSSD.h"
Int_t ver = ITS->IsVersion();
cerr<<"ITS version "<<ver<<" has been found !\n";
- ITS->MakeTreeC();
+ ITS->MakeTreeC();
// Set the models for cluster finding
AliITSgeom *geom = ITS->GetITSgeom();
TClonesArray *recp0 = ITS->ClustersAddress(0);
AliITSClusterFinderSPD *rec0=new AliITSClusterFinderSPD(seg0,dig0,recp0);
ITS->SetReconstructionModel(0,rec0);
+ // test
+ printf("SPD dimensions %f %f \n",seg0->Dx(),seg0->Dz());
+ printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx());
+
// SDD
AliITSDetType *iDetType=ITS->DetType(1);
// end test
// SDD
- // Set response parameters
+ //Set response functions
+ Float_t baseline = 10.;
+ Float_t noise = 1.75;
+
// SDD compression param: 2 fDecrease, 2fTmin, 2fTmax or disable, 2 fTolerance
AliITSDetType *iDetType=ITS->DetType(1);
AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
res1=new AliITSresponseSDD();
ITS->SetResponseModel(1,res1);
}
- Float_t baseline;
- Float_t noise;
- res1->GetNoiseParam(noise,baseline);
- Float_t noise_after_el = res1->GetNoiseAfterElectronics();
- cout << "noise_after_el: " << noise_after_el << endl;
- Float_t fCutAmp;
- fCutAmp = baseline;
- fCutAmp += (2.*noise_after_el); // noise
- cout << "Cut amplitude: " << fCutAmp << endl;
- Int_t cp[8]={0,0,fCutAmp,fCutAmp,0,0,0,0};
- res1->SetCompressParam(cp);
-
- res1->Print();
+ res1->SetMagicValue(900.);
+
+ Float_t maxadc = res1->MaxAdc();
+ Float_t topValue = res1->MagicValue();
+ Float_t norm = maxadc/topValue;
+
+ Float_t fCutAmp = baseline + 2.*noise;
+ fCutAmp *= norm;
+ Int_t cp[8]={0,0,fCutAmp,fCutAmp,0,0,0,0}; //1D
+
+ //res1->SetZeroSupp("2D");
+ res1->SetZeroSupp("1D");
+ res1->SetNoiseParam(noise,baseline);
+ res1->SetDo10to8(kTRUE);
+ res1->SetCompressParam(cp);
+ res1->SetMinVal(4);
+ res1->SetDiffCoeff(3.6,40.);
+ //res1->SetMagicValue(96.95);
+
AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
if (!seg1) {
seg1 = new AliITSsegmentationSDD(geom,res1);
ITS->SetSegmentationModel(1,seg1);
}
AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1);
+ sim1->SetDoFFT(1);
+ sim1->SetCheckNoise(kFALSE);
ITS->SetSimulationModel(1,sim1);
// SSD
AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
ITS->SetSimulationModel(2,sim2);
- cerr<<"Digitizing ITS...\n";
- if(!gAlice->TreeD()) gAlice->MakeTree("D");
- printf("TreeD %p\n",gAlice->TreeD());
- //make branch
- ITS->MakeBranch("D");
+ cerr<<"Digitizing ITS...\n";
TStopwatch timer;
timer.Start();
AliITSresponseSPD::AliITSresponseSPD()
{
// constructor
- SetDiffCoeff();
+ SetThresholds();
SetNoiseParam();
SetDataType();
- SetMinVal();
-
}
//
// Configuration methods
//
- virtual void SetDiffCoeff(Float_t p1=0.00433,Float_t dummy=0.) {
- // Diffusion coefficient
+
+
+ virtual void SetDiffCoeff(Float_t p1=0) {
+ //
fDiffCoeff=p1;
}
- virtual void DiffCoeff(Float_t &diffc,Float_t &dummy) {
- // Get diffusion coefficient
- diffc= fDiffCoeff;
+ virtual Float_t DiffCoeff() {
+ //
+ return fDiffCoeff;
}
- virtual void SetNoiseParam(Float_t n=200., Float_t b=0.) {
- // set noise and baseline
- fNoise=n; fBaseline=b;
- }
- virtual void GetNoiseParam(Float_t &n, Float_t &b) {
- // get noise and baseline
- n=fNoise; b=fBaseline;
- }
- virtual void SetMinVal(Int_t p1=2000) {
- // Zero-suppression option threshold
- fThreshold=p1;
+ virtual void SetThresholds(Float_t thresh=2000, Float_t sigma=280) {
+ // Set Threshold and noise + threshold fluctuations parameter values
+ fThresh=thresh; fSigma=sigma;
}
- virtual Int_t MinVal() {
- // Get zero-suppression threshold
- return fThreshold;
+ virtual void Thresholds(Float_t &thresh, Float_t &sigma) {
+ // Get Threshold and noise + threshold fluctuations parameter values
+ thresh=fThresh; sigma=fSigma;
}
- virtual void SetDataType(const char *data="simulated") {
+ virtual void SetNoiseParam(Float_t col=0., Float_t row=0.) {
+ // set coupling parameters
+ fCouplCol=col; fCouplRow=row;
+ }
+ virtual void GetNoiseParam(Float_t &col, Float_t &row) {
+ // get coupling parameters
+ col=fCouplCol; row=fCouplRow;
+ }
+ virtual void SetDataType(char *data="simulated") {
// Type of data - real or simulated
fDataType=data;
}
- virtual const char *DataType() const {
+ virtual const char *DataType() {
// Get data typer
return fDataType.Data();
}
protected:
- Float_t fDiffCoeff; // Diffusion Coefficient
- Float_t fNoise; // Noise value
- Float_t fBaseline; // Baseline value
- Int_t fThreshold; // Zero-Suppression threshold
-
+ Float_t fDiffCoeff; // Sigma diffusion coefficient (not used)
+ Float_t fThresh; // Threshold value
+ Float_t fSigma; // Noise + threshold fluctuations value
+ Float_t fCouplCol; // Coupling probability along a column
+ Float_t fCouplRow; // Coupling probability along a row
+
TString fDataType; // Type of data - real or simulated
};
#include "AliITShit.h"
#include "AliITSdigit.h"
#include "AliITSmodule.h"
-#include "AliITSMapA2.h"
+#include "AliITSMapA2.h"
#include "AliITSsimulationSPD.h"
#include "AliITSsegmentation.h"
#include "AliITSresponse.h"
-
-
ClassImp(AliITSsimulationSPD)
////////////////////////////////////////////////////////////////////////
// Version: 0
-// Written by Boris Batyunya
-// December 20 1999
+// Written by Rocco Caliandro
+// from a model developed with T. Virgili and R.A. Fini
+// June 15 2000
//
// AliITSsimulationSPD is the simulation of SPDs
+//
//________________________________________________________________________
-
-AliITSsimulationSPD::AliITSsimulationSPD()
-{
+AliITSsimulationSPD::AliITSsimulationSPD(){
// constructor
- fResponse = 0;
+ fResponse = 0;
fSegmentation = 0;
- fMapA2=0;
- fHis = 0;
- fNoise=0.;
- fBaseline=0.;
- fNPixelsZ=0;
- fNPixelsX=0;
+ fHis = 0;
+ fThresh = 0.;
+ fSigma = 0.;
+ fCouplCol = 0.;
+ fCouplRow = 0.;
}
-
-
//_____________________________________________________________________________
AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) {
- // standard constructor
-
- fHis = 0;
+ // constructor
fResponse = resp;
fSegmentation = seg;
- fResponse->GetNoiseParam(fNoise,fBaseline);
-
+ fResponse->Thresholds(fThresh,fSigma);
+ fResponse->GetNoiseParam(fCouplCol,fCouplRow);
+
fMapA2 = new AliITSMapA2(fSegmentation);
-
- //
-
+
+ //
fNPixelsZ=fSegmentation->Npz();
fNPixelsX=fSegmentation->Npx();
-
+ fHis=0;
}
//_____________________________________________________________________________
}
}
-
//__________________________________________________________________________
AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
// Copy Constructor
if(&source == this) return;
this->fMapA2 = source.fMapA2;
- this->fNoise = source.fNoise;
- this->fBaseline = source.fBaseline;
+ this->fThresh = source.fThresh;
+ this->fSigma = source.fSigma;
+ this->fCouplCol = source.fCouplCol;
+ this->fCouplRow = source.fCouplRow;
this->fNPixelsX = source.fNPixelsX;
this->fNPixelsZ = source.fNPixelsZ;
this->fHis = source.fHis;
// Assignment operator
if(&source == this) return *this;
this->fMapA2 = source.fMapA2;
- this->fNoise = source.fNoise;
- this->fBaseline = source.fBaseline;
+ this->fThresh = source.fThresh;
+ this->fSigma = source.fSigma;
+ this->fCouplCol = source.fCouplCol;
+ this->fCouplRow = source.fCouplRow;
this->fNPixelsX = source.fNPixelsX;
this->fNPixelsZ = source.fNPixelsZ;
this->fHis = source.fHis;
}
//_____________________________________________________________________________
-void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy)
-{
+void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module,
+ Int_t dummy) {
// digitize module
- const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
- // for 3.6 eV/pair
- const Float_t kconv = 10000.; // cm -> microns
- Float_t spdLength = fSegmentation->Dz();
- Float_t spdWidth = fSegmentation->Dx();
+ TObjArray *fHits = mod->GetHits();
+ Int_t nhits = fHits->GetEntriesFast();
+ if (!nhits) return;
- Float_t difCoef, dum;
- fResponse->DiffCoeff(difCoef,dum);
- Float_t zPix0 = 1e+6;
- Float_t xPix0 = 1e+6;
- Float_t yPrev = 1e+6;
+ //printf("Module %d (%d hits) \n",module+1,nhits);
- Float_t zPitch = fSegmentation->Dpz(0);
- Float_t xPitch = fSegmentation->Dpx(0);
-
- TObjArray *fHits = mod->GetHits();
- Int_t nhits = fHits->GetEntriesFast();
- if (!nhits) return;
-
- //cout<<"len,wid,dy,nx,nz,pitchx,pitchz ="<<spdLength<<","<<spdWidth<<","<<fSegmentation->Dy()<<","<<fNPixelsX<<","<<fNPixelsZ<<","<<xPitch<<","<<zPitch<<endl;
- // Array of pointers to the label-signal list
-
- Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;;
- Float_t **pList = new Float_t* [maxNDigits];
- memset(pList,0,sizeof(Float_t*)*maxNDigits);
- Int_t indexRange[4] = {0,0,0,0};
-
- // Fill detector maps with GEANT hits
- // loop over hits in the module
- static Bool_t first;
- Int_t lasttrack=-2;
- Int_t hit, iZi, jz, jx;
- //cout<<"SPD: module,nhits ="<<module<<","<<nhits<<endl;
- Int_t idhit=-1;
- for (hit=0;hit<nhits;hit++) {
- AliITShit *iHit = (AliITShit*) fHits->At(hit);
- Int_t layer = iHit->GetLayer();
- Float_t yPix0 = -73;
- if(layer == 1) yPix0 = -77;
-
- if(iHit->StatusEntering()) idhit=hit;
- Int_t itrack = iHit->GetTrack();
- Int_t dray = 0;
-
- if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
-
- // Int_t parent = iHit->GetParticle()->GetFirstMother();
- Int_t partcode = iHit->GetParticle()->GetPdgCode();
-
-// partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
-// 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
-
- /*
- Float_t px = iHit->GetPXL(); // the momenta at the
- Float_t py = iHit->GetPYL(); // each GEANT step
- Float_t pz = iHit->GetPZL();
- Float_t ptot = 1000*sqrt(px*px+py*py+pz*pz);
- */
-
- Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
- // vertex
- pmod *= 1000;
-
-
- if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
- // at p < 6 MeV/c
-
-
- // Get hit z and x(r*phi) cordinates for each module (detector)
- // in local system.
-
- Float_t zPix = kconv*iHit->GetZL();
- Float_t xPix = kconv*iHit->GetXL();
- Float_t yPix = kconv*iHit->GetYL();
-
- // Get track status
- Int_t status = iHit->GetTrackStatus();
- //cout<<"hit,status,y ="<<hit<<","<<status<<","<<yPix<<endl;
-
- // Check boundaries
- if(zPix > spdLength/2) {
- //cout<<"!!!1 z outside ="<<zPix<<endl;
- zPix = spdLength/2 - 10;
- //cout<<"!!!2 z outside ="<<zPix<<endl;
- }
- if(zPix < 0 && zPix < -spdLength/2) {
- //cout<<"!!!1 z outside ="<<zPix<<endl;
- zPix = -spdLength/2 + 10;
- //cout<<"!!!2 z outside ="<<zPix<<endl;
- }
- if(xPix > spdWidth/2) {
- //cout<<"!!!1 x outside ="<<xPix<<endl;
- xPix = spdWidth/2 - 10;
- //cout<<"!!!2 x outside ="<<xPix<<endl;
- }
- if(xPix < 0 && xPix < -spdWidth/2) {
- //cout<<"!!!1 x outside ="<<xPix<<endl;
- xPix = -spdWidth/2 + 10;
- //cout<<"!!!2 x outside ="<<xPix<<endl;
- }
- Int_t trdown = 0;
-
- // enter Si or after event in Si
- if (status == 66 ) {
- zPix0 = zPix;
- xPix0 = xPix;
- yPrev = yPix;
- }
-
- Float_t depEnergy = iHit->GetIonization();
- // skip if the input point to Si
-
- if(depEnergy <= 0.) continue;
-
- // if track returns to the opposite direction:
- if (yPix < yPrev) {
- trdown = 1;
- }
-
-
- // take into account the holes diffusion inside the Silicon
- // the straight line between the entrance and exit points in Si is
- // divided into the several steps; the diffusion is considered
- // for each end point of step and charge
- // is distributed between the pixels through the diffusion.
-
-
- // ---------- the diffusion in Z (beam) direction -------
-
- Float_t charge = depEnergy*kEnToEl; // charge in e-
- Float_t drPath = 0.;
- Float_t tang = 0.;
- Float_t sigmaDif = 0.;
- Float_t zdif = zPix - zPix0;
- Float_t xdif = xPix - xPix0;
- Float_t ydif = TMath::Abs(yPix - yPrev);
- Float_t ydif0 = TMath::Abs(yPrev - yPix0);
-
- if(ydif < 1) continue; // ydif is not zero
-
- Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
-
- Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
- Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
-
- // number of the steps along the track:
- Int_t nsteps = ndZ;
- if(ndX > ndZ) nsteps = ndX;
- if(nsteps < 6) nsteps = 6; // minimum number of the steps
-
- if (projDif < 5 ) {
- drPath = (yPix-yPix0)*1.e-4;
- drPath = TMath::Abs(drPath); // drift path in cm
- sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm
- sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
- nsteps = 1;
- }
-
- if(projDif > 5) tang = ydif/projDif;
- Float_t dCharge = charge/nsteps; // charge in e- for one step
- Float_t dZ = zdif/nsteps;
- Float_t dX = xdif/nsteps;
-
- for (iZi = 1;iZi <= nsteps;iZi++) {
- Float_t dZn = iZi*dZ;
- Float_t dXn = iZi*dX;
- Float_t zPixn = zPix0 + dZn;
- Float_t xPixn = xPix0 + dXn;
-
- if(projDif >= 5) {
- Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
- drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
- if(trdown == 0) {
- drPath = TMath::Abs(drPath) + ydif0*1.e-4;
- }
- if(trdown == 1) {
- drPath = ydif0*1.e-4 - TMath::Abs(drPath);
- drPath = TMath::Abs(drPath);
- }
- sigmaDif = difCoef*sqrt(drPath);
- sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
- }
-
- zPixn = (zPixn + spdLength/2.);
- xPixn = (xPixn + spdWidth/2.);
- Int_t nZpix, nXpix;
- fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
- zPitch = fSegmentation->Dpz(nZpix);
- fSegmentation->GetPadTxz(xPixn,zPixn);
- // set the window for the integration
- Int_t jzmin = 1;
- Int_t jzmax = 3;
- if(nZpix == 1) jzmin =2;
- if(nZpix == fNPixelsZ) jzmax = 2;
-
- Int_t jxmin = 1;
- Int_t jxmax = 3;
- if(nXpix == 1) jxmin =2;
- if(nXpix == fNPixelsX) jxmax = 2;
-
- Float_t zpix = nZpix;
- Float_t dZright = zPitch*(zpix - zPixn);
- Float_t dZleft = zPitch - dZright;
-
- Float_t xpix = nXpix;
- Float_t dXright = xPitch*(xpix - xPixn);
- Float_t dXleft = xPitch - dXright;
-
- Float_t dZprev = 0.;
- Float_t dZnext = 0.;
- Float_t dXprev = 0.;
- Float_t dXnext = 0.;
-
- for(jz=jzmin; jz <=jzmax; jz++) {
- if(jz == 1) {
- dZprev = -zPitch - dZleft;
- dZnext = -dZleft;
- }
- if(jz == 2) {
- dZprev = -dZleft;
- dZnext = dZright;
- }
- if(jz == 3) {
- dZprev = dZright;
- dZnext = dZright + zPitch;
- }
- // kz changes from 1 to the fNofPixels(270)
- Int_t kz = nZpix + jz -2;
-
- Float_t zArg1 = dZprev/sigmaDif;
- Float_t zArg2 = dZnext/sigmaDif;
- Float_t zProb1 = TMath::Erfc(zArg1);
- Float_t zProb2 = TMath::Erfc(zArg2);
- Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
-
-
- // ----------- holes diffusion in X(r*phi) direction --------
-
- if(dZCharge > 1.) {
- for(jx=jxmin; jx <=jxmax; jx++) {
- if(jx == 1) {
- dXprev = -xPitch - dXleft;
- dXnext = -dXleft;
- }
- if(jx == 2) {
- dXprev = -dXleft;
- dXnext = dXright;
- }
- if(jx == 3) {
- dXprev = dXright;
- dXnext = dXright + xPitch;
- }
- Int_t kx = nXpix + jx -2;
-
- Float_t xArg1 = dXprev/sigmaDif;
- Float_t xArg2 = dXnext/sigmaDif;
- Float_t xProb1 = TMath::Erfc(xArg1);
- Float_t xProb2 = TMath::Erfc(xArg2);
- Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
-
- if(dXCharge > 1.) {
- Int_t index = kz-1;
-
- if (first) {
- indexRange[0]=indexRange[1]=index;
- indexRange[2]=indexRange[3]=kx-1;
- first=kFALSE;
- }
-
- indexRange[0]=TMath::Min(indexRange[0],kz-1);
- indexRange[1]=TMath::Max(indexRange[1],kz-1);
- indexRange[2]=TMath::Min(indexRange[2],kx-1);
- indexRange[3]=TMath::Max(indexRange[3],kx-1);
-
- // build the list of digits for this module
- Double_t signal=fMapA2->GetSignal(index,kx-1);
- signal+=dXCharge;
- fMapA2->SetHit(index,kx-1,(double)signal);
- } // dXCharge > 1 e-
- } // jx loop
- } // dZCharge > 1 e-
- } // jz loop
- } // iZi loop
-
- if (status == 65) { // the step is inside of Si
- zPix0 = zPix;
- xPix0 = xPix;
- }
- yPrev = yPix;
- if(dray == 0) {
- GetList(itrack,idhit,pList,indexRange);
- }
+ Int_t number=10000;
+ Int_t *frowpixel = new Int_t[number];
+ Int_t *fcolpixel = new Int_t[number];
+ Double_t *fenepixel = new Double_t[number];
- lasttrack=itrack;
- } // hit loop inside the module
+ // Array of pointers to store the track index of the digits
+ // leave +1, otherwise pList crashes when col=256, row=192
+ Int_t maxNdigits = fNPixelsX*fNPixelsZ+fNPixelsX+1;
+ Float_t **pList = new Float_t* [maxNdigits];
+ memset(pList,0,sizeof(Float_t*)*maxNdigits);
+
+
+ // noise setting
+ SetFluctuations(pList);
-
- // introduce the electronics effects and do zero-suppression
- ChargeToSignal(pList);
- // clean memory
- fMapA2->ClearMap();
+ // loop over hits in the module
+ Int_t hitpos;
+ for (hitpos=0;hitpos<nhits;hitpos++) {
+ HitToDigit(mod,hitpos,module,frowpixel,fcolpixel,fenepixel,pList);
+ }// end loop over digits
+
+ CreateDigit(nhits,module,pList);
+
+ // clean memory
+ delete[] frowpixel;
+ delete[] fcolpixel;
+ delete[] fenepixel;
+ fMapA2->ClearMap();
+ delete [] pList;
+}
+//_____________________________________________________________________________
-}
+void AliITSsimulationSPD::UpdateMap( Int_t row, Int_t col, Double_t ene) {
+//
+// updates the Map of signal, adding the energy (ene) released by the current track
+//
+ Double_t signal;
+ signal = fMapA2->GetSignal(row,col);
+ signal += ene;
+ fMapA2->SetHit(row,col,signal);
+
+ }
+//_____________________________________________________________________________
+void AliITSsimulationSPD::HitToDigit(AliITSmodule *mod, Int_t hitpos, Int_t module,
+ Int_t *frowpixel, Int_t *fcolpixel,
+ Double_t *fenepixel, Float_t **pList) {
+//
+// Steering function to determine the digits associated to a given hit (hitpos)
+// The digits are created by charge sharing (ChargeSharing) and by
+// capacitive coupling (SetCoupling). At all the created digits is associated
+// the track number of the hit (ntrack)
+//
-//---------------------------------------------
-void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
-{
- // lop over nonzero digits
+ static Float_t x1l,y1l,z1l;
+ Float_t x2l,y2l,z2l,etot;
+ Int_t layer,r1,r2,c1,c2,row,col,npixel = 0;
+ Int_t ntrack,idhit;
+ Double_t ene;
+ const Float_t kconv = 10000.; // cm -> microns
+ const Float_t kconv1= 0.277e9; // GeV -> electrons equivalent
+
+ TObjArray *fHits = mod->GetHits();
+ AliITShit *hit = (AliITShit*) fHits->At(hitpos);
+ layer = hit->GetLayer();
+ etot=kconv1*hit->GetIonization();
+ ntrack=hit->GetTrack();
+ idhit=mod->GetHitHitIndex(hitpos);
+
+
+ /*
+ printf("\n layer,etot,ntrack,status %d %f %d %d\n",layer,etot,ntrack,hit->GetTrackStatus()); //debug
+ Int_t idtrack; //debug
+ mod->GetHitTrackAndHitIndex(hitpos,idtrack,idhit);
+ printf("idtrack,idhit %d %d\n",idtrack,idhit); //debug
+ printf("(Dx, Dz)=(%f, %f)\n",fSegmentation->Dx(),fSegmentation->Dz()); //debug
+ */
+
- //set protection
- for(int k=0;k<4;k++) {
- if (indexRange[k] < 0) indexRange[k]=0;
- }
- for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
- for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
+ if (hit->GetTrackStatus()==66) {
+ hit->GetPositionL(x1l,y1l,z1l);
+ // positions shifted and converted in microns
+ x1l = x1l*kconv + fSegmentation->Dx()/2.;
+ z1l = z1l*kconv + fSegmentation->Dz()/2.;
+ //printf("(x1l, z2l)=(%f, %f)\n",x1l,z1l); //debug
+ }
+ else {
+ hit->GetPositionL(x2l,y2l,z2l);
+ // positions shifted and converted in microns
+ x2l = x2l*kconv + fSegmentation->Dx()/2.;
+ z2l = z2l*kconv + fSegmentation->Dz()/2.;
+ //printf("(x2l, z2l)=(%f, %f)\n",x2l,z2l); //debug
+
+
+
+ // to account for the effective sensitive area
+ // introduced in geometry
+ if (z1l<0 || z1l>fSegmentation->Dz()) return;
+ if (z2l<0 || z2l>fSegmentation->Dz()) return;
+ if (x1l<0 || x1l>fSegmentation->Dx()) return;
+ if (x2l<0 || x2l>fSegmentation->Dx()) return;
+
+ //Get the col and row number starting from 1
+ // the x direction is not inverted for the second layer!!!
+ fSegmentation->GetPadIxz(x1l, z1l, c1, r1);
+ fSegmentation->GetPadIxz(x2l, z2l, c2, r2);
+
+ //printf("(c1, r1)=(%d, %d) (c2, r2)=(%d, %d)\n",c1,r1,c2,r2); //debug
+
+ // to account for unexpected equal entrance and
+ // exit coordinates
+ if (x1l==x2l) x2l=x2l+x2l*0.000001;
+ if (z1l==z2l) z2l=z2l+z2l*0.000001;
+
+
+ if ((r1==r2) && (c1==c2))
+ {
+ // no charge sharing
+ npixel = 1;
+ frowpixel[npixel-1] = r1;
+ fcolpixel[npixel-1] = c1;
+ fenepixel[npixel-1] = etot;
+ }
+ else {
+ // charge sharing
+ ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
+ npixel,frowpixel,fcolpixel,fenepixel);
+
+ }
+
+
+ for (Int_t npix=0;npix<npixel;npix++)
+ {
+ row = frowpixel[npix];
+ col = fcolpixel[npix];
+ ene = fenepixel[npix];
+ UpdateMap(row,col,ene);
+ GetList(ntrack,idhit,pList,row,col);
+ // Starting capacitive coupling effect
+ SetCoupling(row,col,ntrack,idhit,pList);
+ }
+ x1l=x2l;
+ y1l=y2l;
+ z1l=z2l;
+ }
+}
- Float_t signal=fMapA2->GetSignal(iz,ix);
+//_________________________________________________________________________
- if (!signal) continue;
+void AliITSsimulationSPD::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
+ Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
+ Int_t r2,Float_t etot,
+ Int_t &npixel,Int_t *frowpixel,
+ Int_t *fcolpixel,Double_t *fenepixel){
+ //
+ // Take into account the geometrical charge sharing when the track
+ // crosses more than one pixel.
+ //
+ //Begin_Html
+ /*
+ <img src="picts/ITS/barimodel_2.gif">
+ </pre>
+ <br clear=left>
+ <font size=+2 color=red>
+ <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
+ </font>
+ <pre>
+ */
+ //End_Html
+
+
+ Float_t xa,za,xb,zb,dx,dz,dtot,dm,refr,refm,refc;
+ Float_t refn=0.;
+ Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
+ Int_t dirx,dirz,rb,cb;
+
+
+ Int_t flag,flagrow,flagcol;
+
+ Double_t epar;
- Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
- if(!pList[globalIndex]){
- //
- // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
- //
+ npixel = 0;
+ xa = x1l;
+ za = z1l;
+ dx = TMath::Abs(x1l-x2l);
+ dz = TMath::Abs(z1l-z2l);
+ dtot = TMath::Sqrt((dx*dx)+(dz*dz));
+ dm = (x2l - x1l) / (z2l - z1l);
- pList[globalIndex] = new Float_t [9];
+ dirx = (Int_t) ((x2l - x1l) / dx);
+ dirz = (Int_t) ((z2l - z1l) / dz);
+
+
+ // calculate the x coordinate of the pixel in the next column
+ // and the z coordinate of the pixel in the next row
- // set list to -3
+ Float_t xpos, zpos;
- *pList[globalIndex] = -3.;
- *(pList[globalIndex]+1) = -3.;
- *(pList[globalIndex]+2) = -3.;
- *(pList[globalIndex]+3) = 0.;
- *(pList[globalIndex]+4) = 0.;
- *(pList[globalIndex]+5) = 0.;
- *(pList[globalIndex]+6) = -1.;
- *(pList[globalIndex]+7) = -1.;
- *(pList[globalIndex]+8) = -1.;
+ fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos);
+ Float_t xsize = fSegmentation->Dpx(0);
+ Float_t zsize = fSegmentation->Dpz(r1-1);
- *pList[globalIndex] = (float)label;
- *(pList[globalIndex]+3) = signal;
- *(pList[globalIndex]+6) = (float)idhit;
- }
- else{
+ if (dirx == 1) refr = xpos+xsize/2.;
+ else refr = xpos-xsize/2.;
- // check the signal magnitude
+ if (dirz == 1) refn = zpos+zsize/2.;
+ else refn = zpos-zsize/2.;
- Float_t highest = *(pList[globalIndex]+3);
- Float_t middle = *(pList[globalIndex]+4);
- Float_t lowest = *(pList[globalIndex]+5);
+
+ flag = 0;
+ flagrow = 0;
+ flagcol = 0;
+ do
+ {
+
+ // calculate the x coordinate of the intersection with the pixel
+ // in the next cell in row direction
+
+ refm = (refn - z1l)*dm + x1l;
+
+ // calculate the z coordinate of the intersection with the pixel
+ // in the next cell in column direction
+
+ refc = (refr - x1l)/dm + z1l;
+
+
+ arefm = refm * dirx;
+ arefr = refr * dirx;
+ arefn = refn * dirz;
+ arefc = refc * dirz;
+
+
+ if ((arefm < arefr) && (arefn < arefc)){
+
+ // the track goes in the pixel in the next cell in row direction
+ xb = refm;
+ zb = refn;
+ cb = c1;
+ rb = r1 + dirz;
+ azb = zb * dirz;
+ az2l = z2l * dirz;
+ if (rb == r2) flagrow=1;
+ if (azb > az2l) {
+ zb = z2l;
+ xb = x2l;
+ }
+
+ // shift to the pixel in the next cell in row direction
+ Float_t zsizeNext = fSegmentation->Dpz(rb-1);
+ //to account for cell at the borders of the detector
+ if(zsizeNext==0) zsizeNext = zsize;
+
+ refn += zsizeNext*dirz;
- signal -= (highest+middle+lowest);
+ }
+ else {
+
+ // the track goes in the pixel in the next cell in column direction
+ xb = refr;
+ zb = refc;
+ cb = c1 + dirx;
+ rb = r1;
+ axb = xb * dirx;
+ ax2l = x2l * dirx;
+ if (cb == c2) flagcol=1;
+ if (axb > ax2l) {
+ zb = z2l;
+ xb = x2l;
+ }
+
+ // shift to the pixel in the next cell in column direction
+ Float_t xsizeNext = fSegmentation->Dpx(cb-1);
+ //to account for cell at the borders of the detector
+ if(xsizeNext==0) xsizeNext = xsize;
+
+ refr += xsizeNext*dirx;
+
+ }
+
+ //calculate the energy lost in the crossed pixel
+ epar = TMath::Sqrt((xb-xa)*(xb-xa)+(zb-za)*(zb-za));
+ epar = etot*(epar/dtot);
+
+ //store row, column and energy lost in the crossed pixel
+ frowpixel[npixel] = r1;
+ fcolpixel[npixel] = c1;
+ fenepixel[npixel] = epar;
+ npixel++;
+
+ // the exit point of the track is reached
+ if (epar == 0) flag = 1;
+ if ((r1 == r2) && (c1 == c2)) flag = 1;
+ if (flag!=1) {
+ r1 = rb;
+ c1 = cb;
+ xa = xb;
+ za = zb;
+ }
+
+ } while (flag==0);
- //
- // compare the new signal with already existing list
- //
+}
+//___________________________________________________________________________
+void AliITSsimulationSPD::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
+ Int_t idhit, Float_t **pList) {
+ //
+ // Take into account the coupling between adiacent pixels.
+ // The parameters probcol and probrow are the fractions of the
+ // signal in one pixel shared in the two adjacent pixels along
+ // the column and row direction, respectively.
+ //
+ //Begin_Html
+ /*
+ <img src="picts/ITS/barimodel_3.gif">
+ </pre>
+ <br clear=left>
+ <font size=+2 color=red>
+ <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
+ </font>
+ <pre>
+ */
+ //End_Html
+
+
+ Int_t j1,j2,flag=0;
+ Double_t pulse1,pulse2;
+
+
+ j1 = row;
+ j2 = col;
+
+ pulse1 = fMapA2->GetSignal(row,col);
+ pulse2 = pulse1;
+
+ for (Int_t isign=-1;isign<=1;isign+=2)
+ {
+
+// loop in row direction
+
+ do
+ {
+ j1 += isign;
+ pulse1 *= fCouplRow;
+
+ if ((j1 < 0) || (j1 > fNPixelsZ-1) || (pulse1 < fThresh))
+ {
+ pulse1 = fMapA2->GetSignal(row,col);
+ j1 = row;
+ flag = 1;
+ }
+ else{
+ UpdateMap(j1,col,pulse1);
+ GetList(ntrack,idhit,pList,j1,col);
+ flag = 0;
+ }
+
+ } while(flag == 0);
+
+
+// loop in column direction
+
+ do
+ {
+ j2 += isign;
+ pulse2 *= fCouplCol;
+
+ if ((j2 < 0) || (j2 > (fNPixelsX-1)) || (pulse2 < fThresh))
+ {
+ pulse2 = fMapA2->GetSignal(row,col);
+ j2 = col;
+ flag = 1;
+ }
+ else{
+ UpdateMap(row,j2,pulse2);
+ GetList(ntrack,idhit,pList,row,j2);
+ flag = 0;
+ }
+
+ } while(flag == 0);
+
+ }
- if(signal<lowest) continue; // neglect this track
+}
+//___________________________________________________________________________
+void AliITSsimulationSPD::CreateDigit(Int_t nhits, Int_t module, Float_t
+**pList) {
+ //
+ // The pixels are fired if the energy deposited inside them is above
+ // the threshold parameter ethr. Fired pixed are interpreted as digits
+ // and stored in the file digitfilename.
+ //
+
+ AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
+
+
+ Int_t digits[3];
+ Int_t tracks[3];
+ Int_t hits[3];
+ Float_t charges[3];
+ Int_t gi,j1;
+
+ if (nhits > 0) {
+
+ for (Int_t r=1;r<=fNPixelsZ;r++) {
+ for (Int_t c=1;c<=fNPixelsX;c++) {
+
+ // check if the deposited energy in a pixel is above the threshold
+ Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
+ gi =r*fNPixelsX+c; // global index
+ if ( signal > fThresh) {
+ digits[0] = r-1; // digits starts from 0
+ digits[1] = c-1; // digits starts from 0
+ //digits[2] = 1;
+ digits[2] = (Int_t) signal; // the signal is stored in electrons
+ for(j1=0;j1<3;j1++){
+ tracks[j1] = (Int_t)(*(pList[gi]+j1));
+ hits[j1] = (Int_t)(*(pList[gi]+j1+6));
+ charges[j1] = 0;
+ }
+ /* debug
+ printf("digits %d %d %d\n",digits[0],digits[1],digits[2]); //debug
+ printf("tracks %d %d %d\n",tracks[0],tracks[1],tracks[2]); //debug
+ printf("hits %d %d %d\n",hits[0],hits[1],hits[2]); //debug
+ */
+ Float_t phys = 0;
+ aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
+ }//endif of threshold condition
+ if(pList[gi]) delete [] pList[gi];
+ }
+ }// enddo on pixels
+ }
+
+}
+//_____________________________________________________________________________
- if (signal>highest){
- *(pList[globalIndex]+5) = middle;
- *(pList[globalIndex]+4) = highest;
- *(pList[globalIndex]+3) = signal;
+void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit, Float_t **pList,
+ Int_t row, Int_t col) {
+ // loop over nonzero digits
- *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
- *(pList[globalIndex]+1) = *pList[globalIndex];
- *pList[globalIndex] = label;
+ Int_t ix = col;
+ Int_t iz = row;
+ Int_t globalIndex;
+ Float_t signal;
+ Float_t highest,middle,lowest;
- *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
- *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
- *(pList[globalIndex]+6) = idhit;
- }
- else if (signal>middle){
- *(pList[globalIndex]+5) = middle;
- *(pList[globalIndex]+4) = signal;
+
+ signal=fMapA2->GetSignal(iz,ix);
- *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
- *(pList[globalIndex]+1) = label;
- *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
- *(pList[globalIndex]+7) = idhit;
- }
- else{
- *(pList[globalIndex]+5) = signal;
- *(pList[globalIndex]+2) = label;
- *(pList[globalIndex]+8) = idhit;
- }
- }
- } // end of loop pixels in x
- } // end of loop over pixels in z
+ globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 1
-}
+ if(!pList[globalIndex])
+ {
+ //
+ // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
+ //
+ pList[globalIndex] = new Float_t [9];
-//---------------------------------------------
-void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
-{
- // add noise and electronics, perform the zero suppression and add the
- // digit to the list
- AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
-
+ // set list to -3
+ *(pList[globalIndex]) = -3.;
+ *(pList[globalIndex]+1) = -3.;
+ *(pList[globalIndex]+2) = -3.;
+ *(pList[globalIndex]+3) = 0.;
+ *(pList[globalIndex]+4) = 0.;
+ *(pList[globalIndex]+5) = 0.;
+ *(pList[globalIndex]+6) = -1.;
+ *(pList[globalIndex]+7) = -1.;
+ *(pList[globalIndex]+8) = -1.;
- Float_t threshold = (float)fResponse->MinVal();
-
- Int_t digits[3], tracks[3], hits[3],gi,j1;
- Float_t charges[3];
- Float_t electronics;
- Float_t signal,phys;
- for(Int_t iz=0;iz<fNPixelsZ;iz++){
- for(Int_t ix=0;ix<fNPixelsX;ix++){
- electronics = fBaseline + fNoise*gRandom->Gaus();
- signal = (float)fMapA2->GetSignal(iz,ix);
- signal += electronics;
- gi =iz*fNPixelsX+ix; // global index
- if (signal > threshold) {
- digits[0]=iz;
- digits[1]=ix;
- digits[2]=1;
- for(j1=0;j1<3;j1++){
- if (pList[gi]) {
- //b.b. tracks[j1]=-3;
- tracks[j1] = (Int_t)(*(pList[gi]+j1));
- hits[j1] = (Int_t)(*(pList[gi]+j1+6));
- }else {
- tracks[j1]=-2; //noise
- hits[j1] = -1;
- }
- charges[j1] = 0;
- }
-
- if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
- tracks[1] = -3;
- hits[1] = -1;
- tracks[2] = -3;
- hits[2] = -1;
- }
- if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
- tracks[1] = -3;
- hits[1] = -1;
- }
- if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
- tracks[2] = -3;
- hits[2] = -1;
- }
- if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
- tracks[2] = -3;
- hits[2] = -1;
- }
-
- phys=0;
- aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
- }
- if(pList[gi]) delete [] pList[gi];
- }
+ *pList[globalIndex] = (float)label;
+ *(pList[globalIndex]+3) = signal;
+ *(pList[globalIndex]+6) = (float)idhit;
}
- delete [] pList;
+ else{
-}
+ // check the signal magnitude
+ highest = *(pList[globalIndex]+3);
+ middle = *(pList[globalIndex]+4);
+ lowest = *(pList[globalIndex]+5);
-//____________________________________________
-void AliITSsimulationSPD::CreateHistograms()
-{
- // create 1D histograms for tests
+ signal -= (highest+middle+lowest);
- printf("SPD - create histograms\n");
+ //
+ // compare the new signal with already existing list
+ //
+ if(signal<lowest) return; // neglect this track
+
+ if (signal>highest)
+ {
+ *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
+ *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
+ *(pList[globalIndex]+6) = idhit;
+ *(pList[globalIndex]+5) = middle;
+ *(pList[globalIndex]+4) = highest;
+ *(pList[globalIndex]+3) = signal;
+ *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
+ *(pList[globalIndex]+1) = *pList[globalIndex];
+ *(pList[globalIndex]) = label;
+ }
+ else if (signal>middle)
+ {
+ *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
+ *(pList[globalIndex]+7) = idhit;
+ *(pList[globalIndex]+5) = middle;
+ *(pList[globalIndex]+4) = signal;
+ *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
+ *(pList[globalIndex]+1) = label;
+ }
+ else
+ {
+ *(pList[globalIndex]+8) = idhit;
+ *(pList[globalIndex]+5) = signal;
+ *(pList[globalIndex]+2) = label;
+ }
+ }
+}
+//_________________________________________________________________________
+void AliITSsimulationSPD::SetFluctuations(Float_t **pList) {
+ //
+ // Set the electronic noise and threshold non-uniformities to all the
+ // pixels in a detector.
+ // The parameter fSigma is the squared sum of the sigma due to noise
+ // and the sigma of the threshold distribution among pixels.
+ //
+ //Begin_Html
+ /*
+ <img src="picts/ITS/barimodel_1.gif">
+ </pre>
+ <br clear=left>
+ <font size=+2 color=red>
+ <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
+ </font>
+ <pre>
+ */
+ //End_Html
+
+
+ Double_t signal;
+
+ Int_t iz,ix;
+ for(iz=1;iz<=fNPixelsZ;iz++){
+ for(ix=1;ix<=fNPixelsX;ix++){
+ signal = fSigma*gRandom->Gaus();
+ fMapA2->SetHit(iz,ix,signal);
+
+ // insert in the label-signal-hit list the pixels fired only by noise
+ if ( signal > fThresh) {
+ Int_t globalIndex = iz*fNPixelsX+ix;
+ pList[globalIndex] = new Float_t [9];
+ *(pList[globalIndex]) = -2.;
+ *(pList[globalIndex]+1) = -2.;
+ *(pList[globalIndex]+2) = -2.;
+ *(pList[globalIndex]+3) = signal;
+ *(pList[globalIndex]+4) = 0.;
+ *(pList[globalIndex]+5) = 0.;
+ *(pList[globalIndex]+6) = -1.;
+ *(pList[globalIndex]+7) = -1.;
+ *(pList[globalIndex]+8) = -1.;
+ }
+ } // end of loop on pixels
+ } // end of loop on pixels
+
+ }
+//____________________________________________
+
+void AliITSsimulationSPD::CreateHistograms() {
+ // CreateHistograms
+
+ Int_t i;
fHis=new TObjArray(fNPixelsZ);
- for (Int_t i=0;i<fNPixelsZ;i++) {
- TString spdName("spd_");
- Char_t pixelz[4];
- sprintf(pixelz,"%d",i+1);
- spdName.Append(pixelz);
- (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
+ for(i=0;i<fNPixelsZ;i++) {
+ TString spdname("spd_");
+ Char_t candnum[4];
+ sprintf(candnum,"%d",i+1);
+ spdname.Append(candnum);
+ (*fHis)[i] = new TH1F(spdname.Data(),"SPD maps",
fNPixelsX,0.,(Float_t) fNPixelsX);
}
+
}
//____________________________________________
-void AliITSsimulationSPD::ResetHistograms()
-{
+void AliITSsimulationSPD::ResetHistograms() {
//
// Reset histograms for this detector
//
-
- for ( int i=0;i<fNPixelsZ;i++ ) {
+ Int_t i;
+ for(i=0;i<fNPixelsZ;i++ ) {
if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
}
}
-
-
-
-
-
-
-
-
-
AliITSsimulationSPD(const AliITSsimulationSPD &source); // copy constructor
AliITSsimulationSPD& operator=(const AliITSsimulationSPD &source); // ass. operator
- void DigitiseModule(AliITSmodule *mod,Int_t module,Int_t dummy);
- void ChargeToSignal(Float_t **pList);
- void GetList(Int_t track, Int_t hit, Float_t **pList, Int_t *IndexRange);
+ void DigitiseModule(AliITSmodule *mod,Int_t module, Int_t dummy);
+ void SetFluctuations(Float_t **pList);
+ void HitToDigit(AliITSmodule *mod, Int_t hitpos, Int_t module,
+ Int_t *frowpixel, Int_t *fcolpixel, Double_t *fenepixel,
+ Float_t **pList);
+
+ void UpdateMap( Int_t row, Int_t col, Double_t ene);
+ void ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
+ Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
+ Int_t r2,Float_t etot,
+ Int_t &npixel,Int_t *frowpixel,
+ Int_t *fcolpixel,Double_t *fenepixel);
+
+ void SetCoupling(Int_t row, Int_t col, Int_t ntrack, Int_t idhit, Float_t **pList);
+ void CreateDigit(Int_t nhits, Int_t module, Float_t **pList);
+ void GetList(Int_t track,Int_t idhit, Float_t **pList, Int_t row, Int_t col);
void CreateHistograms();
void ResetHistograms();
private:
AliITSMapA2 *fMapA2; // MapA2
- Float_t fNoise; // Noise
- Float_t fBaseline; // Baseline
+ Float_t fThresh; // Threshold
+ Float_t fSigma; // Noise
+ Float_t fCouplCol; // Coupling along columns
+ Float_t fCouplRow; // Coupling along rows
Int_t fNPixelsX; // NPixelsX
Int_t fNPixelsZ; // NPixelsZ
#endif
-
-
//Test ITS reconstruction
gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClusters.C");
+
+ delete gAlice; gAlice=0;
+
if (rc=AliITSFindClusters()) return rc;
//gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSgraphycs.C");
// simulation but in cluster finder as well, please set them via your
// local Config.C - the streamer will take care of writing the correct
// info and you'll no longer be obliged to set them again for your cluster
- // finder as it's done in this macro
+ // finder as it's done in this macro (ugly and impractical, no? )
// SPD
-
-
ITS->MakeTreeC();
Int_t nparticles=gAlice->GetEvent(0);
+
AliITSDetType *iDetType=ITS->DetType(0);
AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
TClonesArray *dig0 = ITS->DigitsAddress(0);
TClonesArray *recp0 = ITS->ClustersAddress(0);
AliITSClusterFinderSPD *rec0=new AliITSClusterFinderSPD(seg0,dig0,recp0);
ITS->SetReconstructionModel(0,rec0);
+ // test
+ //printf("SPD dimensions %f %f \n",seg0->Dx(),seg0->Dz());
+ //printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx());
+
// SDD
AliITSDetType *iDetType=ITS->DetType(2);
AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
- seg2->SetDetSize(72960.,40000.,303.);
TClonesArray *dig2 = ITS->DigitsAddress(2);
AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
ITS->SetReconstructionModel(2,rec2);
for (int nev=evNumber1; nev<= evNumber2; nev++) {
if(nev>0) {
- nparticles = gAlice->GetEvent(nev);
- gAlice->SetEvent(nev);
- if(!gAlice->TreeR()) gAlice-> MakeTree("R");
- ITS->MakeBranch("R");
+ nparticles = gAlice->GetEvent(nev);
+ gAlice->SetEvent(nev);
+ if(!gAlice->TreeR()) gAlice-> MakeTree("R");
+ ITS->MakeBranch("R");
}
cout << "nev " <<nev<<endl;
cout << "nparticles " <<nparticles<<endl;
Int_t last_entry=0;
timer.Start();
ITS->DigitsToRecPoints(nev,last_entry,"All");
+ //ITS->DigitsToRecPoints(nev,last_entry,"SPD");
timer.Stop(); timer.Print();
} // event loop
file->Close();
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
ITS->SetResponseModel(1,res1);
}
+
//res1->SetChargeLoss(0.);
Float_t baseline;
Float_t noise;
res1->GetNoiseParam(noise,baseline);
Float_t noise_after_el = res1->GetNoiseAfterElectronics();
- cout << "noise_after_el: " << noise_after_el << endl;
+ cout << "noise_after_el: " << noise_after_el << endl;
Float_t fCutAmp;
fCutAmp = baseline;
fCutAmp += (2.*noise_after_el); // noise
AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1);
ITS->SetSimulationModel(1,sim1);
sim1->Print();
-
-
+
// SPD
AliITSDetType *iDetType=ITS->DetType(0);
AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
AliITSresponseSPD *res0 = (AliITSresponseSPD*)iDetType->GetResponseModel();
+ //AliITSresponseSPD *res0= new AliITSresponseSPD();
+ //ITS->SetResponseModel(0,res0);
+
+// to change the parameters
+ //res0->SetThresholds(2000,280);
+ //res0->SetNoiseParam(0., 0.);
+ //res0->SetNoiseParam(0.04, 0.08);
+
+// to monitor the parameters
+ Float_t thresh, sigma;
+ res0->Thresholds(thresh, sigma);
+ printf("SPD: threshold %e sigma %e\n",thresh, sigma);
+ Float_t col, row;
+ res0->GetNoiseParam(col, row);
+ printf("SPD: Coupling by column %e Coupling by row %e\n",col, row);
+
AliITSsimulationSPD *sim0=new AliITSsimulationSPD(seg0,res0);
ITS->SetSimulationModel(0,sim0);
- // test
- //printf("SPD dimensions %f %f \n",seg0->Dx(),seg0->Dz());
- //printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx());
- //printf("SPD pitches %d %d \n",seg0->Dpz(0),seg0->Dpx(0));
- // end test
-
// SSD
AliITSDetType *iDetType=ITS->DetType(2);
AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
- seg2->SetDetSize(72960.,40000.,303.);
AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
res2->SetSigmaSpread(3.,2.);
AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
for (Int_t nev=evNumber1; nev<= evNumber2; nev++) {
cout << "nev " <<nev<<endl;
if(nev>0) {
- nparticles = gAlice->GetEvent(nev);
- gAlice->SetEvent(nev);
- if(!gAlice->TreeD()) gAlice-> MakeTree("D");
- ITS->MakeBranch("D");
+ nparticles = gAlice->GetEvent(nev);
+ gAlice->SetEvent(nev);
+ if(!gAlice->TreeD()) gAlice-> MakeTree("D");
+ ITS->MakeBranch("D");
}
cout << "nparticles " <<nparticles<<endl;
if (nev < evNumber1) continue;
if(nsignal) nbgr_ev=Int_t(nev/nsignal);
timer.Start();
ITS->HitsToDigits(nev,nbgr_ev,size," ","All"," ");
+ //ITS->HitsToDigits(nev,nbgr_ev,size," ","SPD"," ");
timer.Stop(); timer.Print();
} // event loop
- delete sim0;
+// delete sim0;
delete sim1;
delete sim2;
-
file->Close();
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
#pragma link C++ class AliITSsegmentationSSD+;
#pragma link C++ class AliITSresponse+;
#pragma link C++ class AliITSresponseSPD+;
-#pragma link C++ class AliITSresponseSPDbari+;
+#pragma link C++ class AliITSresponseSPDdubna+;
#pragma link C++ class AliITSresponseSDD+;
#pragma link C++ class AliITSresponseSSD+;
#pragma link C++ class AliITSsimulation+;
#pragma link C++ class AliITSsimulationSPD+;
-#pragma link C++ class AliITSsimulationSPDbari+;
+#pragma link C++ class AliITSsimulationSPDdubna+;
#pragma link C++ class AliITSsimulationSDD+;
#pragma link C++ class AliITSsimulationSSD+;
#pragma link C++ class AliITSsimulationFastPoints+;
#pragma link C++ class AliITSsimulationFastPointsV0+;
#pragma link C++ class AliITSClusterFinder+;
#pragma link C++ class AliITSClusterFinderSPD+;
-#pragma link C++ class AliITSClusterFinderSPDbari+;
+#pragma link C++ class AliITSClusterFinderSPDdubna+;
#pragma link C++ class AliITSClusterFinderSDD+;
#pragma link C++ class AliITSClusterFinderSSD+;
#pragma link C++ class AliITSDetType+;
AliITSGeant3Geometry.cxx \
AliITSsimulationFastPoints.cxx \
AliITSsimulationFastPointsV0.cxx AliITSsimulation.cxx \
- AliITSsimulationSPD.cxx AliITSsimulationSPDbari.cxx \
+ AliITSsimulationSPD.cxx AliITSsimulationSPDdubna.cxx \
AliITSsimulationSDD.cxx \
AliITSetfSDD.cxx AliITSsimulationSSD.cxx AliITSdcsSSD.cxx \
AliITSdigit.cxx AliITSRawCluster.cxx AliITSRecPoint.cxx \
AliITSsegmentation.cxx AliITSsegmentationSPD.cxx \
AliITSsegmentationSDD.cxx AliITSsegmentationSSD.cxx\
AliITSresponse.cxx AliITSresponseSPD.cxx \
- AliITSresponseSPDbari.cxx \
+ AliITSresponseSPDdubna.cxx \
AliITSresponseSDD.cxx AliITSresponseSSD.cxx \
AliITSClusterFinder.cxx AliITSClusterFinderSPD.cxx \
- AliITSClusterFinderSPDbari.cxx \
+ AliITSClusterFinderSPDdubna.cxx \
AliITSClusterFinderSDD.cxx AliITSRawData.cxx \
AliITSHuffman.cxx AliITSClusterFinderSSD.cxx \
AliITSclusterSSD.cxx AliITSpackageSSD.cxx \
-#include "iostream.h"
-
-void SPDclusterTest (Int_t evNumber1=0,Int_t evNumber2=0)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- } else {
- delete gAlice;
- gAlice=0;
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!file) file = new TFile("galice.root");
- file->ls();
-
-// Get AliRun object from file or create it if not on file
+void SPDclusterTest(Int_t evNumber1=0,Int_t evNumber2=0){
+ //
+ // macro to monitor the SPD digitization and clusterization done with
+ // the Bari/Salerno model
+ //
+ // R. Caliandro 15/05/2001
+ //
+ //
+ //--plots displayed:
+ //
+ //--pag1: number of hits per SPD detector (1-->250)
+ // number of hits per SPD detector (1-->250)
+ // number of clusters per SPD detector (1-->250)
+ //
+ //--pag2: r-phi cluster length layer 1 (red)
+ // z cluster length layer 1 (red)
+ // r-phi cluster length layer 2 (blue)
+ // z cluster length layer 2 (blue)
+ //
+ //--pag3: r-phi resolution layer 1 (red)
+ // z resolution layer 1 (red)
+ // r-phi resolution layer 2 (blue)
+ // z resolution layer 2 (blue)
+ //
+ //--pag4: Cluster shape analysis for clusters of 1, 2 and 3 digits
+ // zdim versus xdim for clusters of 4 digits
+ //
+ // input file name, digitized and clusterized
+ char *filein="galice.root";
+ // output file name, containing histograms
+ char *fileout="SPD_his.root";
+ // flag for debugging: 0=no debugging, 1=debugging
+ Int_t debug=0;
+
+
+ // Dynamically link some shared libs
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } else {
+ delete gAlice;
+ gAlice=0;
+ }
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+
+ // Connect the Root Galice file containing Geometry, Kine and Hits
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filein);
+ if (!file) file = new TFile(filein);
+
+ // Get AliRun object from file or create it if not on file
+ if (!gAlice) {
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (gAlice) printf("AliRun object found on file\n");
+ if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+ }
+
+ //to get the segmentation pointer
+ AliITS *ITS = (AliITS*) gAlice->GetModule("ITS");
+ AliITSDetType *iDetType=ITS->DetType(0);
+ AliITSsegmentationSPD *seg=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
+
+//=======================================================
+//--booking of ntuples
+
+//--ntuple for each detector
+ TNtuple *ntuple2 = new TNtuple("ntuple2","","ndet:lay:lad:det:nhits:ndig:nclus");
+//--ntuple for each cluster
+ TNtuple *ntuple = new TNtuple("ntuple","","ndet:iclus:ndigclus:xdim:zdim:xdiff:zdiff:anglex:anglez:pmom:errx:errz");
+
+//--booking of histograms
+//layer 1
+ TH1F *hist1n1 = new TH1F("hist1n1","xdim",15,0.5,15.5);
+ TH1F *hist2n1 = new TH1F("hist2n1","zdim",10,0.5,10.5);
+ TH1F *hist3n1 = new TH1F("hist3n1","dig/clus",20,0.5,20.5);
+ TH1F *hist4n1 = new TH1F("hist4n1","errx",100,0,0.01);
+ TH1F *hist5n1 = new TH1F("hist5n1","errz",500,0,0.05);
+ TH2F *hist7n1 = new TH2F("hist7n1","xdim:delx",80,0,800.,15,0.5,15.5);
+ TH2F *hist8n1 = new TH2F("hist8n1","zdim:delz",180,0,1800.,10,0.5,10.5);
+//layer 2
+ TH1F *hist1n2 = new TH1F("hist1n2","xdim",15,0.5,15.5);
+ TH1F *hist2n2 = new TH1F("hist2n2","zdim",10,0.5,10.5);
+ TH1F *hist3n2 = new TH1F("hist3n2","dig/clus",20,0.5,20.5);
+ TH1F *hist4n2 = new TH1F("hist4n2","errx",100,0,0.01);
+ TH1F *hist5n2 = new TH1F("hist5n2","errz",500,0,0.05);
+ TH2F *hist7n2 = new TH2F("hist7n2","xdim:delx",80,0,800.,15,0.5,15.5);
+ TH2F *hist8n2 = new TH2F("hist8n2","zdim:delz",180,0,1800.,10,0.5,10.5);
+//--resolution
+ TH1F *hist1 = new TH1F("hist1","xdiff",200,-100,100);
+ TH1F *hist3 = new TH1F("hist3","xdiff",200,-100,100);
+ TH1F *hist2 = new TH1F("hist2","zdiff",170,-850,850);
+ TH1F *hist4 = new TH1F("hist4","zdiff",170,-850,850);
+//--momentum
+ TH1F *hist5 = new TH1F("hist5","pmom",200,0,2000);
+//--rapidity
+ TH1F *hist6 = new TH1F("hist6","rapidity",60,-3,3);
+ TH1F *hist6b= new TH1F("hist6b","rapidity - charged tracks",60,-3,3);
+ TH1F *hist6b1= new TH1F("hist6b1","rapidity - charged tracks SPD",60,-3,3);
+//--pseudo-rapidity
+ TH1F *hist6p = new TH1F("hist6p","eta - charged tracks ",60,-3,3);
+ TH1F *hist6p1 = new TH1F("hist6p1","eta - charged tracks SPD ",60,-3,3);
+ TH1F *hist6p2 = new TH1F("hist6p2","eta - charged tracks SPD 2 ",60,-3,3);
+//--resolution vs angle
+ TH1F *hist11n1=new TH1F("hist11n1","anglex - layer 1",180,-90,90);
+ TH1F *hist11n2=new TH1F("hist11n2","anglex - layer 2",180,-90,90);
+ TH1F *hist12n1=new TH1F("hist12n1","anglez - layer 1",360,-180,180);
+ TH1F *hist12n2=new TH1F("hist12n2","anglez - layer 2",360,-180,180);
+ TH2F *hist13n1=new TH2F("hist13n1","xidff:anglex",20,-15,15,200,-100,100);
+ TH2F *hist13n2=new TH2F("hist13n2","xidff:anglex",20,-30,-15,200,-100,100);
+ TH2F *hist14n1=new TH2F("hist14n1","zidff:anglez",18,-90,90,170,-850,850);
+ TH2F *hist14n2=new TH2F("hist14n2","zidff:anglez",18,-90,90,170,-850,850);
+//--histograms for cluster shape analysis
+ TH1F *histsp1=new TH1F("histsp1","Cluster shape (1)",10,0.5,10.5);
+ TH2F *histsp2=new TH2F("histsp2","Cluster shape (2)",5,0.5,5.5,5,0.5,5.5);
+//=======================================================
+
+ //loop over events
+ for (int nev=0; nev<= evNumber2; nev++) {
+ Int_t nparticles = gAlice->GetEvent(nev);
+ cout << "nev " <<nev<<endl;
+ cout << "nparticles " <<nparticles<<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks = TH->GetEntries();
+ cout << "ntracks " <<ntracks<<endl;
+
+ // Get pointers to Alice detectors and Digit containers
+ AliITS *ITS = (AliITS *)gAlice->GetModule("ITS");
+ TClonesArray *Particles = gAlice->Particles();
+ if(!ITS) return;
+
+ // fill modules with sorted by module hits
+ Int_t nmodules;
+ ITS->InitModules(-1,nmodules);
+ ITS->FillModules(nev,evNumber2,nmodules," "," ");
+
+ // get pointer to modules array
+ TObjArray *mods = ITS->GetModules();
+ AliITShit *itsHit;
+
+
+ //get the Tree for clusters
+ ITS->GetTreeC(nev);
+ TTree *TC=ITS->TreeC();
+ //TC->Print();
+ Int_t nent=TC->GetEntries();
+ printf("Found %d entries in the tree of clusters)\n",nent);
+ TClonesArray *ITSclusters = ITS->ClustersAddress(0);
+ printf("ITSclusters %p\n",ITSclusters);
+
+ //get the Tree for digits
+ TTree *TD = gAlice->TreeD();
+ //TD->Print();
+ Int_t nentd=TD->GetEntries();
+ printf("Found %d entries in the tree of digits)\n",nentd);
+ TObjArray *fBranches=TD->GetListOfBranches();
+ TBranch *branch = (TBranch*)fBranches->UncheckedAt(0);
+ printf ("branch %p entries %d \n",branch,branch->GetEntries());
+ TClonesArray *ITSdigits = ITS->DigitsAddress(0);
+ printf ("ITSdigits %p \n",ITSdigits);
+
+ //get the Tree for rec points
+ TTree *TR = gAlice->TreeR();
+ //TR->Print();
+ Int_t nentr=TR->GetEntries();
+ printf("Found %d entries in the tree of rec points)\n",nentr);
+ TClonesArray *ITSrec = ITS->RecPoints();
+ printf ("ITSrec %p \n",ITSrec);
+ AliITSRecPoint *recp;
+
+ // calculus of rapidity distribution for the generated tracks
+ gAlice-> ResetHits();
+ TParticle *particle;
+ for (Int_t track=0; track<ntracks; track++)
+ {
+ particle = (TParticle*)gAlice->Particle(track);
+ Int_t ikparen = particle -> GetFirstMother();
+ Double_t charge = particle -> GetPDG() ->Charge();
+ charge = charge/3.; //charge is multiplied by 3 in PDG
+ Double_t mass = particle -> GetPDG() -> Mass();
+ Double_t eta = particle -> Eta();
+ Int_t pdgcode = particle -> GetPdgCode();
+ char* title = particle -> GetTitle();
+ if (ikparen<0)
+ {
+ Double_t part_ene = particle->Energy();
+ Double_t part_pz = particle->Pz();
+ Double_t rapid;
+ if (part_ene != part_pz)
+ {
+ rapid=0.5*TMath::Log((part_ene+part_pz)/(part_ene-part_pz));
+ }
+ else {
+ rapid = 1.e30;
+ }
+ // filling of the rapidity histogram
+ hist6->Fill( (Float_t) rapid);
+ if( charge != 0 ) {
+ hist6b->Fill( (Float_t) rapid);
+ hist6p->Fill( (Float_t) eta);
+ }
+// printf("charge= %f, mass = %f , pdg= %d, title = %s\n",
+// charge,mass,pdgcode,title);
+ }
+ }
+
+ AliITSgeom *g = ((AliITS *)ITS)->GetITSgeom();
+ Int_t lay, lad, det;
+ //printf("Starts loop on SPD detectors\n");
+
+
+ //loop over the pixel detectors index=0-79 (1-20)*4 layer 1
+ // index=80-239 (1-40)*4 layer 2
+ for (Int_t index=g->GetStartSPD();index<=g->GetLastSPD();index++)
+// for (Int_t index=g->GetStartSPD();index<1;index++) //debug
+ {
+
+ g->GetModuleId(index,lay,lad,det);
+ //printf("detector %d (lay=%d lad=%d det=%d)\n",index+1,lay,lad,det);
+
+ AliITSmodule *itsModule = (AliITSmodule*) mods->At(index);
+ Int_t numofhits = itsModule->GetNhits();
+ //printf("number of hits %d\n",numofhits);
+ if(!numofhits) continue;
+
+ //---------- starts test on digits
+ ITS->ResetDigits();
+ TD->GetEvent(index);
+ Int_t ndigits = ITSdigits->GetEntriesFast();
+ //if (ndigits) printf("Found %d digits for module %d \n",ndigits,index+1);
+ if (!ndigits) printf("no digits found \n");
+
+
+
+ if(debug==1) {
+ //loop on digits
+ for (Int_t digit=0;digit<ndigits;digit++) {
+ ITSdigit = (AliITSdigitSPD*)ITSdigits->UncheckedAt(digit);
+ printf("digit=%d fCoord1=%d FCoord2=%d fSignal=%d fTracks=%d fHits=%d \n",digit,ITSdigit->fCoord1,ITSdigit->fCoord2,ITSdigit->fSignal,ITSdigit->fTracks[0],ITSdigit->fHits[0]);
+ }
+ cout<<"END test for digits "<<endl;
+ }
+
+
+ //---------- starts test on clusters
+ ITS->ResetClusters();
+ TC->GetEvent(index);
+ Int_t nclust = ITSclusters->GetEntries();
+ //printf("Found %d clusters \n",nclust);
+ if (!nclust) printf("no clusters found \n");
+
+
+ if(debug==1) {
+ //loop on clusters
+ for (Int_t clu=0;clu<nclust;clu++)
+ {
+ //itsclu = (AliITSRawClusterSPD*) ITSclusters->UncheckedAt(clu);
+ itsclu = (AliITSRawClusterSPD*) ITSclusters->At(clu);
+ printf("cluster %d nZ=%f nX=%f Q=%f Z=%f X=%f\n",clu+1,itsclu->NclZ(),
+ itsclu->NclX(),itsclu->Q(),itsclu->Z(),itsclu->X());
+ }
+ cout<<"END test for clusters "<<endl;
+ }
+
+
+
+ //---------- starts test on rec points
+ ITS->ResetRecPoints();
+ TR->GetEvent(index);
+ Int_t nrecpoints = ITSrec->GetEntries();
+ //printf("Found %d recpoints for module %d \n",nrecpoints,index+1);
+ if (!nrecpoints) printf("no recpoints found \n");
+
+ if(debug==1) {
+ //loop on rec points
+ for (Int_t irec=0;irec<nrecpoints;irec++) {
+ recp = (AliITSRecPoint*)ITSrec->UncheckedAt(irec);
+ printf("%d %f %f %f %f %d %d %d\n",irec+1,recp->GetX(),recp->GetZ(),
+ recp->fSigmaX2,recp->fSigmaZ2,
+ recp->fTracks[0],recp->fTracks[1],recp->fTracks[2]);
+ }
}
- // ------------ Cluster and point analysis histogramms ------------
-
- TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.);
- TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.);
- TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.);
- TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.);
- TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.);
- TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.);
- TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.);
- TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.);
-
- TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.);
- TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.);
- TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.);
- TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.);
-
- TH1F *Ptot1 = new TH1F("Ptot1","Total momentum (GeV/C) for layers 1",100,0.,5.);
- TH1F *Pz1 = new TH1F("Pz1","Pz (GeV/C) for layers 1",100,-5.,5.);
- TH1F *Theta1 = new TH1F("Theta1","Theta angle (rad) for layers 1",100,0.,4.);
- TH1F *Y1 = new TH1F("Y1","Rapidity for layers 1",100,-4.,4.);
- TH1F *Eta1 = new TH1F("Eta1","PseudoRapidity for layers 1",100,-4.,4.);
- TH1F *Y1Den = new TH1F("Y1Den","Rapidity for layers 1",100,-0.5,0.5);
- TH1F *Eta1Den = new TH1F("Eta1Den","PseudoRapidity for layers 1",100,-0.5,0.5);
- TH1F *Y1DenA = new TH1F("Y1DenA","Rapidity for layers 1",100,-0.5,0.5);
- TH1F *Eta1DenA = new TH1F("Eta1DenA","PseudoRapidity for layers 1",100,-0.5,0.5);
- TH1F *Phi1 = new TH1F("Phi1","Phi angle (rad) for layers 1",100,0.,7.);
- TH1F *Ptot2 = new TH1F("Ptot2","Total momentum (GeV/C) for layers 2",100,0.,5.);
- TH1F *Pz2 = new TH1F("Pz2","Pz (GeV/C) for layers 2",100,-5.,5.);
- TH1F *Theta2 = new TH1F("Theta2","Theta angle (rad) for layers 2",100,0.,4.);
- TH1F *Y2 = new TH1F("Y2","Rapidity for layers 2",100,-4.,4.);
- TH1F *Eta2 = new TH1F("Eta2","PseudoRapidity for layers 2",100,-4.,4.);
- TH1F *Y2Den = new TH1F("Y2Den","Rapidity for layers 2",100,-0.5,0.5);
- TH1F *Eta2Den = new TH1F("Eta2Den","PseudoRapidity for layers 2",100,-0.5,0.5);
- TH1F *Y2DenA = new TH1F("Y2DenA","Rapidity for layers 2",100,-0.5,0.5);
- TH1F *Eta2DenA = new TH1F("Eta2DenA","PseudoRapidity for layers 2",100,-0.5,0.5);
- TH1F *Phi2 = new TH1F("Phi2","Phi angle (rad) for layers 2",100,0.,7.);
-
- // -------------- Create ntuples --------------------
- // ntuple structures:
-
- struct {
- Int_t lay;
- Int_t nx;
- Int_t nz;
- Int_t hitprim;
- Int_t partcode;
- Int_t ntrover;
- Float_t dx;
- Float_t dz;
- Float_t pmod;
- } ntuple_st;
-
- struct {
- Int_t lay;
- Int_t lad;
- Int_t det;
- Int_t nx;
- Int_t nz;
- Int_t ntrover;
- Int_t noverlaps;
- Int_t noverprim;
- Float_t qcl;
- Float_t dx;
- Float_t dz;
- Float_t x;
- Float_t z;
- } ntuple1_st;
-
- struct {
- // Int_t lay;
- Int_t lay;
- Int_t nx;
- Int_t nz;
- Float_t x;
- Float_t z;
- Float_t qcl;
- } ntuple2_st;
-
- ntuple = new TTree("ntuple","Demo ntuple");
- ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
- ntuple->Branch("nx",&ntuple_st.nx,"nx/I");
- ntuple->Branch("nz",&ntuple_st.nz,"nz/I");
- ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
- ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
- ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
- ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
- ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
- ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
-
- ntuple1 = new TTree("ntuple1","Demo ntuple1");
- ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
- ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
- ntuple1->Branch("det",&ntuple1_st.det,"det/I");
- ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I");
- ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I");
- ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F");
- ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
- ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
- ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
- ntuple1->Branch("x",&ntuple1_st.x,"x/F");
- ntuple1->Branch("z",&ntuple1_st.z,"z/F");
- ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
- ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
-
-
- ntuple2 = new TTree("ntuple2","Demo ntuple2");
- // ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
- ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
- ntuple2->Branch("x",&ntuple2_st.x,"x/F");
- ntuple2->Branch("z",&ntuple2_st.z,"z/F");
- ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I");
- ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I");
- ntuple2->Branch("qcl",&ntuple2_st.qcl,"qcl/F");
-
-// ------------------------------------------------------------------------
+ printf("Detector n.%d (%d hits) (%d digits) (%d clusters)\n",
+ index+1,numofhits,ndigits,nclust);
+
+ // fill ntuple2
+ ntuple2->Fill ( (Float_t) index+1,
+ (Float_t) lay,
+ (Float_t) lad,
+ (Float_t) det,
+ (Float_t) numofhits,
+ (Float_t) ndigits,
+ (Float_t) nclust);
+
+ Int_t xlow;
+ Int_t zlow;
+ Int_t xhigh;
+ Int_t zhigh;
+ Int_t colcenter;
+ Int_t rowcenter;
+
+ // loop on clusters in each detector
+ for (Int_t i=0; i<nclust; i++)
+ {
+
+ irawclu = (AliITSRawClusterSPD*) ITSclusters->UncheckedAt(i);
+ irecp = (AliITSRecPoint*)ITSrec->UncheckedAt(i);
+
+ Int_t xdim = irawclu->NclX();
+ Int_t zdim = irawclu->NclZ();
+ Float_t errx = TMath::Sqrt(irecp->fSigmaX2);
+ Float_t errz = TMath::Sqrt(irecp->fSigmaZ2);
+ Float_t xcenter = irawclu->X();
+ Float_t zcenter = irawclu->Z();
+ Float_t ndigclus = irawclu->Q();
+ Int_t itrackclus = irecp->fTracks[0];
+
+ //Find the hits associated to the main track of the cluster
+ // loop on hits in the detector
+ for (Int_t hit=0; hit<numofhits; hit++)
+ {
+ AliITShit *itsHit = (AliITShit*)itsModule->GetHit(hit);
+ Int_t itrackhit = itsHit->GetTrack();
+ //Take the same track index of the main track of the cluster
+ if (itrackhit == itrackclus) {
+ if (itsHit->GetTrackStatus()==66) {
+ Float_t x1l = 10000*itsHit->GetXL(); //in microns
+ Float_t y1l = 10000*itsHit->GetYL();
+ Float_t z1l = 10000*itsHit->GetZL();
+ Float_t p1x = 1000*itsHit->GetPXL(); //in MeV/c
+ Float_t p1y = 1000*itsHit->GetPYL();
+ Float_t p1z = 1000*itsHit->GetPZL();
+ }
+ else {
+ Float_t x2l = 10000*itsHit->GetXL();
+ Float_t y2l = 10000*itsHit->GetYL();
+ Float_t z2l = 10000*itsHit->GetZL();
+ Float_t p2x = 1000*itsHit->GetPXL();
+ Float_t p2y = 1000*itsHit->GetPYL();
+ Float_t p2z = 1000*itsHit->GetPZL();
+
+
+ }
+ }
+ }// end loop on hits on detector
+
+
+ Float_t pmom=TMath::Sqrt(p1x*p1x+p1y*p1y+p1z*p1z);
+ hist5->Fill(pmom);
+
+ Float_t dxhit = TMath::Abs(x2l-x1l);
+ Float_t dzhit = TMath::Abs(z2l-z1l);
+
+ Float_t xmidhit = (x1l + x2l)/2;
+ Float_t zmidhit = (z1l + z2l)/2;
+
+// printf("cluster n.%d: x=%f z=%f\n",i,xcenter,zcenter);
+// printf("track n.%d: x1=%f x2=%f z1=%f z2=%f\n",itrackclus,
+// x1l, x2l, z1l, z2l);
+
+ // analysis of resolution vs angle
+ if(index<80)
+ {
+ Float_t px = -p1x;
+ Float_t py = -p1y;
+ }
+ else{
+ Float_t px = p1x;
+ Float_t py = p1y;
+ }
+ Float_t pz = p1z;
+ // anglex is the angle in xy plane (local frame)
+ Float_t anglex = atan2(px,py);
+ // anglez is the angle in zy plane (local frame)
+ Float_t anglez = atan2(pz,py);
+ anglex *= 180.0/TMath::Pi(); // degrees
+ anglez *= 180.0/TMath::Pi(); // degrees
+
+ if(xmidhit != 0 || zmidhit != 0)
+ {
+ Float_t xdiff = (xcenter - xmidhit);
+ Float_t zdiff = (zcenter - zmidhit);
+
+ if(index<80)
+ {
+ // resolution plots
+ hist1->Fill(xdiff);
+ hist2->Fill(zdiff);
+
+ // plots of resolution vs angle
+ hist11n1->Fill(anglex);
+ hist12n1->Fill(anglez);
+ hist13n1->Fill(anglex,xdiff);
+ hist14n1->Fill(anglez,zdiff);
+
+ } else {
+
+ // resolution plots
+ hist3->Fill(xdiff);
+ hist4->Fill(zdiff);
+
+ // plots of resolution vs angle
+ hist11n2->Fill(anglex);
+ hist12n2->Fill(anglez);
+ hist13n2->Fill(anglex,xdiff);
+ hist14n2->Fill(anglez,zdiff);
+
+ }
+ }
+
+ // fill the ntuple
+ ntuple->Fill ( (Float_t) index,
+ (Float_t) i,
+ (Float_t) ndigclus,
+ (Float_t) xdim,
+ (Float_t) zdim,
+ (Float_t) xdiff,
+ (Float_t) zdiff,
+ (Float_t) anglex,
+ (Float_t) anglez,
+ (Float_t) pmom,
+ errx,
+ errz);
+
+ // other histograms
+ if(index<80)
+ {
+ hist1n1->Fill((Float_t) xdim);
+ hist2n1->Fill((Float_t) zdim);
+ hist3n1->Fill(ndigclus);
+ hist4n1->Fill(errx);
+ hist5n1->Fill(errz);
+ hist7n1->Fill(dxhit,(Float_t) xdim);
+ hist8n1->Fill(dzhit,(Float_t) zdim);
+
+ } else {
+ hist1n2->Fill((Float_t) xdim);
+ hist2n2->Fill((Float_t) zdim);
+ hist3n2->Fill(ndigclus);
+ hist4n2->Fill(errx);
+ hist5n2->Fill(errz);
+ hist7n2->Fill(dxhit,(Float_t) xdim);
+ hist8n2->Fill(dzhit,(Float_t) zdim);
+ }
+
+ //histograms for cluster shape analysis
+ Int_t xx;
+ if(ndigclus<=3) {
+ if(ndigclus==1) {
+ xx=1;
+ } elseif(ndigclus==2){
+ if(zdim==2 && xdim==1) xx=2;
+ if(zdim==1 && xdim==2) xx=3;
+ if(zdim==2 && xdim==2) xx=4;
+ } elseif(ndigclus==3){
+ if(zdim==2 && xdim==2) xx=5;
+ if(zdim==3 && xdim==1) xx=6;
+ if(zdim==1 && xdim==3) xx=7;
+ if(zdim==3 && xdim==3) xx=8;
+ if(zdim==3 && xdim==2) xx=9;
+ if(zdim==2 && xdim==3) xx=10;
+ }
+ histsp1->Fill((Float_t) xx);
+ } elseif(ndigclus==4){
+ histsp2->Fill((Float_t) xdim,(Float_t) zdim);
+ }
+
+
+ }//end loop on clusters
+
+ } //end loop on the SPD detectors
+
+} //end loop on events
+
+
+//================== Plot the results ===================================
+
+gROOT->Reset();
+gStyle->SetFillColor(0);
+gStyle->SetStatW(0.37);
+gStyle->SetStatH(0.22);
+
+//----------------------------------------------------- page 1
+gStyle->SetOptLogy(0);
+gStyle->SetOptStat(1100);
+
+ TCanvas *c1 = new TCanvas("c1","hits, digits, clusters",200,10,700,780);
+c1->SetFillColor(0);
+c1->Divide(1,3);
+ c1->cd(1);
+ ntuple2->SetMarkerColor(kRed);
+ ntuple2->SetMarkerStyle(20);
+ ntuple2->SetMarkerSize(0.6);
+ ntuple2->Draw("nhits:ndet","");
+ c1->cd(2);
+ ntuple2->SetMarkerColor(kBlue);
+ ntuple2->Draw("ndig:ndet","");
+ c1->cd(3);
+ ntuple2->SetMarkerColor(6);
+ ntuple2->Draw("nclus:ndet","");
+
+//----------------------------------------------------- page 2
+
+ TCanvas *c2 = new TCanvas("c2","Cluster Lengths",200,10,700,780);
+ //
+ // Inside this canvas, we create 4 pads
+ //
+ pad1 = new TPad("pad1","xdim layer1" ,0.01,0.51,0.49,0.99,10);
+ pad2 = new TPad("pad2","zdim layer1" ,0.51,0.51,0.99,0.99,10);
+ pad3 = new TPad("pad3","xdim layer2" ,0.01,0.01,0.49,0.49,10);
+ pad4 = new TPad("pad4","zdim layer2" ,0.51,0.01,0.99,0.49,10);
+ pad1->Draw();
+ pad2->Draw();
+ pad3->Draw();
+ pad4->Draw();
//
-// Loop over events
+ gStyle->SetStatW(0.40);
+ gStyle->SetStatH(0.20);
+ gStyle->SetStatColor(42);
+ gStyle->SetOptStat(111110);
+// gStyle->SetOptFit(1);
+
+ pad1->cd();
+ pad1->SetLogy();
+ hist1n1->Draw();
+ hist1n1->SetLineWidth(2);
+ hist1n1->SetLineColor(kRed);
+ hist1n1->GetXaxis()->SetNdivisions(110);
+ hist1n1->GetYaxis()->SetLabelSize(0.06);
+ c2->Update();
+ //
+ pad2->cd();
+ pad2->SetLogy();
+ hist2n1->Draw();
+ hist2n1->SetLineWidth(2);
+ hist2n1->SetLineColor(kRed);
+ hist2n1->GetXaxis()->SetNdivisions(110);
+ hist2n1->GetYaxis()->SetLabelSize(0.06);
+ c2->Update();
+ //
+ pad3->cd();
+ pad3->SetLogy();
+ hist1n2->Draw();
+ hist1n2->SetLineWidth(2);
+ hist1n2->SetLineColor(kBlue);
+ hist1n2->GetXaxis()->SetNdivisions(110);
+ hist1n2->GetYaxis()->SetLabelSize(0.06);
+ c2->Update();
+ //
+ pad4->cd();
+ pad4->SetLogy();
+ hist2n2->Draw();
+ hist2n2->SetLineColor(kBlue);
+ hist2n2->SetLineWidth(2);
+ hist2n2->GetXaxis()->SetNdivisions(110);
+ hist2n2->GetYaxis()->SetLabelSize(0.06);
+ c2->Update();
+ //
+//----------------------------------------------------- page 3
+
+ TCanvas *c3 = new TCanvas("c3","Resolutions",200,10,700,780);
+ //
+ // Inside this canvas, we create 4 pads
+ //
+ pad1 = new TPad("pad1","xdiff layer1" ,0.01,0.51,0.49,0.99,10);
+ pad2 = new TPad("pad2","zdiff layer1" ,0.51,0.51,0.99,0.99,10);
+ pad3 = new TPad("pad3","xdiff layer2" ,0.01,0.01,0.49,0.49,10);
+ pad4 = new TPad("pad4","zdiff layer2" ,0.51,0.01,0.99,0.49,10);
+ pad1->Draw();
+ pad2->Draw();
+ pad3->Draw();
+ pad4->Draw();
//
- for (int nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nev " << nev <<endl;
- cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
- cout<<"ntracks "<<ntracks<<endl;
-
-// Get pointers to Alice detectors and Digits containers
- AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
- TClonesArray *Particles = gAlice->Particles();
-
- if (ITS) {
- // fill modules with sorted by module hits
- Int_t nmodules;
- ITS->InitModules(-1,nmodules);
- // ITS->FillModules(nev,-1,evNumber2,nmodules," "," ");
- ITS->FillModules(nev,evNumber2,nmodules," "," ");
- //get pointer to modules array
- TObjArray *ITSmodules = ITS->GetModules();
- AliITShit *itsHit;
-
- // get the Tree for clusters
- ITS->GetTreeC(nev);
- TTree *TC=ITS->TreeC();
- Int_t nent=TC->GetEntries();
- printf("Found %d entries in the tree (must be one per module per event!)\n",nent);
- Int_t lay, lad, det;
- AliITSgeom *geom = ITS->GetITSgeom();
-
- for (Int_t idettype=0;idettype<3;idettype++) {
-
- TClonesArray *ITSclusters = ITS->ClustersAddress(idettype);
- //printf ("ITSclusters %p \n",ITSclusters);
-
- if (idettype != 0) continue;
-
- Float_t occup1 = 0;
- Float_t occup2 = 0;
-
- // Module loop
- for (Int_t mod=0; mod<nent; mod++) {
- AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod);
- geom->GetModuleId(mod,lay,lad,det);
-
- Int_t nhits = itsModule->GetNhits();
- //if(nhits) printf("module nhits %d %d\n",mod,nhits);
- if(!nhits) continue;
-
- ITS->ResetClusters();
- TC->GetEvent(mod);
- Int_t nclust = ITSclusters->GetEntries();
- if (!nclust) continue;
-
- // cluster/hit loops
- //cout<<"mod,lay,nclust,nhits ="<<mod<<","<<lay<<","<<nclust<<","<<nhits<<endl;
- for (Int_t clu=0;clu<nclust;clu++) {
- itsclu = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu);
-
- Int_t noverlaps = 0;
- Int_t noverprim = 0;
-
- Int_t clustersizex = itsclu->NclX();
- Int_t clustersizez = itsclu->NclZ();
- // Int_t xstart = itsclu->XStart();
- // Int_t xstop = itsclu->XStop();
- Int_t xstart = itsclu->XStartf();
- Int_t xstop = itsclu->XStopf();
- Float_t fxstart = xstart*50;
- Float_t fxstop = (xstop+1)*50;
- Float_t zstart = itsclu->ZStart();
- Float_t zstop = itsclu->ZStop();
- Int_t zend = itsclu->Zend();
- Int_t ntrover = itsclu->NTracks();
- Float_t clusterx = itsclu->X();
- Float_t clusterz = itsclu->Z();
- Float_t clusterQ = itsclu->Q();
-
- if(lay == 1) occup1 += clusterQ;
- if(lay == 2) occup2 += clusterQ;
-
- ntuple2_st.lay = lay;
- ntuple2_st.x = clusterx/1000.;
- ntuple2_st.z = clusterz/1000.;
- ntuple2_st.nx = clustersizex;
- ntuple2_st.nz = clustersizez;
- ntuple2_st.qcl = clusterQ;
-
- ntuple2->Fill();
-
- Int_t icl = 0;
- Float_t dxprimlast = 10.e+6;
- Float_t dzprimlast = 10.e+6;
-
- Float_t SPDlength = 83600;
- Float_t SPDwidth = 12800;
- Float_t xhit0 = 1e+5;
- Float_t zhit0 = 1e+5;
-
- for (Int_t hit=0;hit<nhits;hit++) {
-
- // Find coordinate differences between the hit and cluster positions
- // for the resolution determination.
-
- itsHit = (AliITShit*)itsModule->GetHit(hit);
-
- Int_t hitlayer = itsHit->GetLayer();
- Int_t hitladder= itsHit->GetLadder();
- Int_t hitdet= itsHit->GetDetector();
- Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
- Int_t track = itsHit->GetTrack();
- Int_t dray = 0;
- Int_t hitstat = itsHit->GetTrackStatus();
-
- Float_t zhit = 10000*itsHit->GetZL();
- Float_t xhit = 10000*itsHit->GetXL();
-
- Float_t pxsimL = itsHit->GetPXL(); // the momenta at GEANT points
- Float_t pysimL = itsHit->GetPYL();
- Float_t pzsimL = itsHit->GetPZL();
- Float_t psimL = TMath::Sqrt(pxsimL*pxsimL+pysimL*pysimL+pzsimL*pzsimL);
-
- // Check boundaries
- if(zhit > SPDlength/2) {
- //cout<<"!!! z outside ="<<zhit<<endl;
- zhit = SPDlength/2 - 10;
- }
- if(zhit < 0 && zhit < -SPDlength/2) {
- //cout<<"!!! z outside ="<<zhit<<endl;
- zhit = -SPDlength/2 + 10;
- }
- if(xhit > SPDwidth/2) {
- //cout<<"!!! x outside ="<<xhit<<endl;
- xhit = SPDwidth/2 - 10;
- }
- if(xhit < 0 && xhit < -SPDwidth/2) {
- //cout<<"!!! x outside ="<<xhit<<endl;
- xhit = -SPDwidth/2 + 10;
- }
-
- zhit += SPDlength/2;
- xhit += SPDwidth/2;
- Float_t yhit = 10000*itsHit->GetYL();
-
- if(hitlayer == 1 && hitstat == 66 && yhit > 71) {
- xhit0 = xhit;
- zhit0 = zhit;
- }
- if(hitlayer == 2 && hitstat == 66 && yhit < -71) {
- xhit0 = xhit;
- zhit0 = zhit;
- }
-
- if(hitstat != 68) continue; // Take only the hit if the last
- // track point went out from
- // the detector.
-
- if(xhit0 > 9e+4 || zhit0 > 9e+4) continue;
-
- Float_t xmed = (xhit + xhit0)/2;
- Float_t zmed = (zhit + zhit0)/2;
-
- Float_t xdif = xmed - clusterx;
- Float_t zdif = zmed - clusterz;
-
-
- // Consider the hits inside of cluster region only
-
- if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
- icl = 1;
-
- // part = (TParticle *)particles.UncheckedAt(track);
- // Int_t partcode = part->GetPdgCode();
- // Int_t primery = gAlice->GetPrimary(track);
-
- Int_t parent = itsHit->GetParticle()->GetFirstMother();
- Int_t partcode = itsHit->GetParticle()->GetPdgCode();
-
-// partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
-// 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
-
- Float_t pmod = itsHit->GetParticle()->P(); // total momentum at the
- // vertex
- Float_t energy = itsHit->GetParticle()->Energy(); // energy at the
- // vertex
- Float_t mass = itsHit->GetParticle()->GetMass(); // particle mass
-
- Float_t pz = itsHit->GetParticle()->Pz(); // z momentum componetnt
- // at the vertex
- Float_t px = itsHit->GetParticle()->Px(); // z momentum componetnt
- // at the vertex
- Float_t py = itsHit->GetParticle()->Py(); // z momentum componetnt
- // at the vertex
- Float_t phi = itsHit->GetParticle()->Phi(); // Phi angle at the
- // vertex
- Float_t theta = itsHit->GetParticle()->Theta(); // Theta angle at the
- // vertex
- //Float_t eta = itsHit->GetParticle()->Eta(); // Pseudo rapidity at the
- // vertex
- if((energy-pz) > 0) {
- Float_t y = 0.5*TMath::Log((energy+pz)/(energy-pz));
- }else{
- cout<<" Warning: energy < pz ="<<energy<<","<<pz<<endl;
- y = 10;
- }
- Float_t eta = -TMath::Log(TMath::Tan(theta/2));
- pmod *= 1.0e+3;
-
-
- Float_t pxsim = itsHit->GetPXG(); // the momenta at this GEANT point
- Float_t pysim = itsHit->GetPYG();
- Float_t pzsim = itsHit->GetPZG();
- Float_t psim = TMath::Sqrt(pxsim*pxsim+pysim*pysim+pzsim*pzsim);
-
- Int_t hitprim = 0;
-
- if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
- // at p < 6 MeV/c
-
- if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but
- // not for delta ray which
- // also went out from the
- // detector and returned
- // again
-
- if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
- // particles
-
- if(hitprim > 0) noverprim = noverprim + 1;
-
- if(hitprim > 0) {
- dxprimlast = xdif;
- dzprimlast = zdif;
- }
-
- // fill ntuple
-
- ntuple_st.lay = hitlayer;
- ntuple_st.nx = clustersizex;
- ntuple_st.nz = clustersizez;
- ntuple_st.hitprim = hitprim;
- ntuple_st.partcode = partcode;
- ntuple_st.ntrover = ntrover;
- ntuple_st.dx = xdif;
- ntuple_st.dz = zdif;
- ntuple_st.pmod = pmod;
-
- ntuple->Fill();
-
- if(hitlayer == 1) {
- Y1DenA->Fill(y);
- Eta1DenA->Fill(eta);
- }
- if(hitlayer == 2) {
- Y2DenA->Fill(y);
- Eta2DenA->Fill(eta);
- }
-
-
-
- if(hitprim > 0) { // for primary particles
-
- if(hitlayer == 1) {
- Xres1->Fill(xdif);
- Zres1->Fill(zdif);
- Ptot1->Fill(pmod/1000.);
- Pz1->Fill(pz);
- Theta1->Fill(theta);
- Y1->Fill(y);
- Eta1->Fill(eta);
- Y1Den->Fill(y);
- Eta1Den->Fill(eta);
- Phi1->Fill(phi);
- }
- if(hitlayer == 2) {
- Xres2->Fill(xdif);
- Zres2->Fill(zdif);
- Ptot2->Fill(pmod/1000.);
- Pz2->Fill(pz);
- Theta2->Fill(theta);
- Y2->Fill(y);
- Eta2->Fill(eta);
- Y2Den->Fill(y);
- Eta2Den->Fill(eta);
- Phi2->Fill(phi);
- }
- } // primery particles
-
- } // end of cluster region
- } // end of hit loop
-
- if(icl == 1) {
-
- // fill ntuple1
-
- // ntuple1->Fill(clusterlayer,clustersizex,clustersizez,noverlaps,\
-noverprim,dx,dz);
-
- if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
- // delta rays only
-
- ntuple1_st.lay = lay;
- ntuple1_st.lad = lad;
- ntuple1_st.det = det;
- ntuple1_st.x = clusterx*1000.;
- ntuple1_st.z = clusterz*1000.;
- ntuple1_st.nx = clustersizex;
- ntuple1_st.nz = clustersizez;
- ntuple1_st.qcl = clusterQ;
- ntuple1_st.ntrover = ntrover;
- ntuple1_st.noverlaps = noverlaps;
- ntuple1_st.noverprim = noverprim;
- ntuple1_st.dx = dxprimlast;
- ntuple1_st.dz = dzprimlast;
-
- ntuple1->Fill();
-
- } // icl = 1
- } // cluster loop
- } // module loop
-
- cout<<" Occupancy for layer-1 ="<<occup1<<endl;
- cout<<" Occupancy for layer-2 ="<<occup2<<endl;
- // The real occupancy values are:
- // (for full ALICE event at the full SPD acceptence)
- // occup1 /= 3932160;
- // occup2 /= 7864320;
-
- } // idettype loop
- } // end if ITS
-} // event loop
-
-
- // Write and Draw Histogramms and ntuples
-
-
-
- TFile fhistos("SPD_his.root","RECREATE");
-
- ntuple->Write();
- ntuple1->Write();
- ntuple2->Write();
-
- Nxpix1->Write();
- Nzpix1->Write();
- Nxpix2->Write();
- Nzpix2->Write();
-
- Xpix1->Write();
- Zpix1->Write();
- Xpix2->Write();
- Zpix2->Write();
-
- Xres1->Write();
- Zres1->Write();
- Xres2->Write();
- Zres2->Write();
-
- Ptot1->Write();
- Pz1->Write();
- Theta1->Write();
- Y1->Write();
- Eta1->Write();
- Y1Den->Write();
- Eta1Den->Write();
- Y1DenA->Write();
- Eta1DenA->Write();
- Phi1->Write();
-
- Ptot2->Write();
- Pz2->Write();
- Theta2->Write();
- Y2->Write();
- Eta2->Write();
- Y2Den->Write();
- Eta2Den->Write();
- Y2DenA->Write();
- Eta2DenA->Write();
- Phi2->Write();
-
- fhistos.Close();
- cout<<"!!! Histogramms and ntuples were written"<<endl;
-
- TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
- c1->Divide(2,2);
- c1->cd(1);
- gPad->SetFillColor(33);
- Xres1->SetFillColor(42);
- Xres1->Draw();
- c1->cd(2);
- gPad->SetFillColor(33);
- Zres1->SetFillColor(46);
- Zres1->Draw();
- c1->cd(3);
- gPad->SetFillColor(33);
- Xres2->SetFillColor(42);
- Xres2->Draw();
- c1->cd(4);
- gPad->SetFillColor(33);
- Zres2->SetFillColor(46);
- Zres2->Draw();
-
- cout<<"END test for clusters and hits "<<endl;
-
-// file->Close();
+ gStyle->SetStatW(0.20);
+ gStyle->SetStatH(0.20);
+ gStyle->SetStatColor(42);
+ gStyle->SetOptStat(0);
+ gStyle->SetOptFit(1);
+
+ pad1->cd();
+ hist1->Draw();
+ hist1->SetLineColor(kRed);
+ hist1->Fit("gaus");
+ c3->Update();
+ //
+ pad2->cd();
+ hist2->Draw();
+ hist2->SetLineColor(kRed);
+ hist2->Fit("gaus");
+ c3->Update();
+ //
+ pad3->cd();
+ hist3->Draw();
+ hist3->SetLineColor(kBlue);
+ hist3->Fit("gaus");
+ c3->Update();
+ //
+ pad4->cd();
+ hist4->Draw();
+ hist4->SetLineColor(kBlue);
+ hist4->Fit("gaus");
+ c3->Update();
+
+//----------------------------------------------------- page 4
+ TCanvas *c4 = new TCanvas("c4","Cluster Shape Analysis",200,10,700,780);
+//
+gStyle->SetOptStat(0);
+c4->SetFillColor(0);
+c4->Divide(1,2);
+
+ c4->cd(1);
+ c4_1->SetLogy();
+ histsp1->Draw();
+ c4->cd(2);
+ c4_2->Divide(2,1);
+ c4_2->cd(1);
+ histsp2->Draw("box");
+ c4_2->cd(2);
+ clustershape();
+
+
+
+
+//================== Store the histograms ===================================
+
+ //to write the histograms into a .root file
+ TFile outfile(fileout,"RECREATE");
+
+ ntuple->Write();
+ ntuple2->Write();
+ hist1n1->Write();
+ hist2n1->Write();
+ hist3n1->Write();
+ hist4n1->Write();
+ hist5n1->Write();
+ hist7n1->Write();
+ hist8n1->Write();
+ hist1n2->Write();
+ hist2n2->Write();
+ hist3n2->Write();
+ hist4n2->Write();
+ hist5n2->Write();
+ hist7n2->Write();
+ hist8n2->Write();
+ hist1->Write();
+ hist2->Write();
+ hist3->Write();
+ hist4->Write();
+ hist5->Write();
+ hist6->Write();
+ hist6b->Write();
+ hist6p->Write();
+ hist6b1->Write();
+ hist6p1->Write();
+ hist6p2->Write();
+ hist11n1->Write();
+ hist11n2->Write();
+ hist12n1->Write();
+ hist12n2->Write();
+ hist13n1->Write();
+ hist13n2->Write();
+ hist14n1->Write();
+ hist14n2->Write();
+ histsp1->Write();
+ histsp2->Write();
+
+ outfile->Close();
+}
+//-----------------------------------------------------------------
+
+
+
+void clustershape(){
+ //
+ //macro to display the legend of the cluster shape analysis plot
+ //
+
+ TPad *pad1 = new TPad("pad1", "This is pad2",0,0,1,1);
+ pad1->Draw();
+ pad1->cd();
+ pad1->Range(0,0.25,1,1);
+ pad1->SetFillColor(0);
+ pad1->SetBorderSize(1);
+
+//------------------------------------------
+ Float_t yfirst= 0.95;
+ Float_t ysh = 0.05;
+ Float_t ysiz = 0.02;
+ Float_t ynum = 0.005;
+//------------------------------------------
+
+ //bin 1
+ TLatex *tex = new TLatex(0.12,yfirst,"1");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 2
+ yfirst=yfirst-ysh;
+ TLatex *tex = new TLatex(0.12,yfirst,"2");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 3
+ yfirst=yfirst-ysh;
+ TLatex *tex = new TLatex(0.12,yfirst-ynum,"3");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst-ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 4
+ yfirst=yfirst-1.8*ysh;
+ TLatex *tex = new TLatex(0.12,yfirst+3*ynum,"4");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 5
+ yfirst=yfirst-1.5*ysh;
+ TLatex *tex = new TLatex(0.12,yfirst+3*ynum,"5");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 6
+ yfirst=yfirst-1.5*ysh;
+ TLatex *tex = new TLatex(0.12,yfirst+ynum,"6");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.7,yfirst,0.9,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 7
+ yfirst=yfirst-2*ysh;
+ TLatex *tex = new TLatex(0.12,yfirst+ysiz+ynum,"7");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.3,yfirst+ysiz,0.5,yfirst+2*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.3,yfirst+2*ysiz,0.5,yfirst+3*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 8
+ yfirst=yfirst-1.5*ysh;
+ TLatex *tex = new TLatex(0.12,yfirst+ynum,"8");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.7,yfirst+2*ysiz,0.9,yfirst+3*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 9
+ yfirst=yfirst-1.5*ysh;
+ TLatex *tex = new TLatex(0.12,yfirst+ynum,"9");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst,0.7,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.7,yfirst+ysiz,0.9,yfirst+2*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+
+ //bin 10
+ yfirst=yfirst-1.5*ysh;
+ TLatex *tex = new TLatex(0.12,yfirst-ynum,"10");
+ tex->SetTextSize(0.07);
+ tex->SetLineWidth(2);
+ tex->Draw();
+ TPave *pave = new TPave(0.3,yfirst-ysiz,0.5,yfirst,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.3,yfirst,0.5,yfirst+ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
+ TPave *pave = new TPave(0.5,yfirst+ysiz,0.7,yfirst+2*ysiz,1,"br");
+ pave->SetFillColor(18);
+ pave->SetLineWidth(2);
+ pave->Draw();
}
-
-
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff