Bari/Salerno model set as defaault SPD simulation
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 31 May 2001 20:37:56 +0000 (20:37 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 31 May 2001 20:37:56 +0000 (20:37 +0000)
13 files changed:
ITS/AliITS.cxx
ITS/AliITSFindClusters.C
ITS/AliITSHits2Digits.C
ITS/AliITSresponseSPD.cxx
ITS/AliITSresponseSPD.h
ITS/AliITSsimulationSPD.cxx
ITS/AliITSsimulationSPD.h
ITS/AliITStest.C
ITS/ITSDigitsToClusters.C
ITS/ITSHitsToDigits.C
ITS/ITSLinkDef.h
ITS/Makefile
ITS/SPDclusterTest.C

index 978ef19..cc8bd9d 100644 (file)
 
 /*
 $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)
 
@@ -231,7 +237,7 @@ the AliITS class.
 #include "AliITSresponse.h"
 #include "AliITSsegmentationSPD.h"
 #include "AliITSresponseSPD.h"
-#include "AliITSresponseSPDbari.h"
+#include "AliITSresponseSPDdubna.h"
 #include "AliITSsegmentationSDD.h"
 #include "AliITSresponseSDD.h"
 #include "AliITSsegmentationSSD.h"
index 7873aae..504fbe3 100644 (file)
@@ -23,7 +23,7 @@ Int_t AliITSFindClusters() {
    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();
 
@@ -34,6 +34,10 @@ Int_t AliITSFindClusters() {
    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);
index 3961f98..2994c05 100644 (file)
@@ -45,7 +45,10 @@ Int_t AliITSHits2Digits()
   // 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();
@@ -53,25 +56,33 @@ Int_t AliITSHits2Digits()
     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
@@ -82,12 +93,8 @@ Int_t AliITSHits2Digits()
   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();
index 0db68b7..18f3d44 100644 (file)
@@ -23,10 +23,8 @@ ClassImp(AliITSresponseSPD)
 AliITSresponseSPD::AliITSresponseSPD()
 {
   // constructor
-   SetDiffCoeff();
+   SetThresholds();
    SetNoiseParam();
    SetDataType();
-   SetMinVal();
-
 }
 
index 1ab6075..dd01821 100644 (file)
@@ -19,35 +19,37 @@ public:
   //
   // 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();
   } 
@@ -56,11 +58,12 @@ public:
     
     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
 };
 
index bec9375..f478d73 100644 (file)
 #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;
 }
 
 //_____________________________________________________________________________
@@ -76,14 +69,15 @@ AliITSsimulationSPD::~AliITSsimulationSPD() {
   }                
 }
 
-
 //__________________________________________________________________________
 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;
@@ -96,8 +90,10 @@ AliITSsimulationSPD&
   //    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;
@@ -105,521 +101,621 @@ AliITSsimulationSPD&
   }
 //_____________________________________________________________________________
 
-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();
     }
 
 }
-
-
-
-
-
-
-
-
-
index 2e022b5..d44f015 100644 (file)
@@ -20,9 +20,22 @@ public:
   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();
@@ -34,8 +47,10 @@ public:
 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
 
@@ -47,5 +62,3 @@ private:
 
 #endif 
 
-
-
index 3351b8c..eb5fe14 100644 (file)
@@ -23,6 +23,9 @@ Int_t AliITStest() {
 
 //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");
index a4937a0..ac61105 100644 (file)
@@ -50,7 +50,7 @@ void ITSDigitsToClusters (Int_t evNumber1=0,Int_t evNumber2=0)
    // 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? )
 
 
 
@@ -58,17 +58,20 @@ void ITSDigitsToClusters (Int_t evNumber1=0,Int_t evNumber2=0)
 
    // 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
 
@@ -99,7 +102,6 @@ void ITSDigitsToClusters (Int_t evNumber1=0,Int_t evNumber2=0)
 
    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);
@@ -123,10 +125,10 @@ void ITSDigitsToClusters (Int_t evNumber1=0,Int_t evNumber2=0)
 
    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;
@@ -136,6 +138,7 @@ void ITSDigitsToClusters (Int_t evNumber1=0,Int_t evNumber2=0)
        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 
 
@@ -145,17 +148,3 @@ void ITSDigitsToClusters (Int_t evNumber1=0,Int_t evNumber2=0)
 
    file->Close();
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
index 4b950bd..b12736a 100644 (file)
@@ -57,12 +57,13 @@ void ITSHitsToDigits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal  =25, In
          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
@@ -86,28 +87,36 @@ void ITSHitsToDigits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal  =25, In
    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);
@@ -136,10 +145,10 @@ void ITSHitsToDigits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal  =25, In
    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;
@@ -149,27 +158,13 @@ void ITSHitsToDigits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal  =25, In
        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();
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
index b60cd31..4dad32e 100644 (file)
 #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+;
index d3f2f59..921f5b6 100644 (file)
@@ -20,7 +20,7 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                 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 \
@@ -28,10 +28,10 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.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 \
index 7076a76..4fd361b 100644 (file)
-#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();
 }
-
-