removed; now the Bari/Salerno model is the default
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Jun 2001 09:50:08 +0000 (09:50 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Jun 2001 09:50:08 +0000 (09:50 +0000)
12 files changed:
ITS/AliITSClusterFinderSPDbari.cxx [deleted file]
ITS/AliITSClusterFinderSPDbari.h [deleted file]
ITS/AliITSFindClustersBari.C [deleted file]
ITS/AliITSHits2DigitsBari.C [deleted file]
ITS/AliITSresponseSPDbari.cxx [deleted file]
ITS/AliITSresponseSPDbari.h [deleted file]
ITS/AliITSsimulationSPDbari.cxx [deleted file]
ITS/AliITSsimulationSPDbari.h [deleted file]
ITS/AliITStestBari.C [deleted file]
ITS/ITSDigitsToClustersBari.C [deleted file]
ITS/ITSHitsToDigitsBari.C [deleted file]
ITS/SPDclusterTestBari.C [deleted file]

diff --git a/ITS/AliITSClusterFinderSPDbari.cxx b/ITS/AliITSClusterFinderSPDbari.cxx
deleted file mode 100644 (file)
index fc1209f..0000000
+++ /dev/null
@@ -1,438 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-
-#include "AliITSClusterFinderSPDbari.h"
-#include "AliITS.h"
-#include "AliITSgeom.h"
-#include "AliITSdigit.h"
-#include "AliITSRawCluster.h"
-#include "AliITSRecPoint.h"
-#include "AliITSsegmentation.h"
-#include "AliITSresponse.h"
-#include "AliRun.h"
-
-
-
-
-ClassImp(AliITSClusterFinderSPDbari)
-
-//----------------------------------------------------------
-AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari
-(AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
-{
-  // constructor
-    fSegmentation=seg;
-    fDigits=digits;
-    fClusters=recp;
-    fNclusters= fClusters->GetEntriesFast();
-    SetDx();
-    SetDz();
-}
-
-//_____________________________________________________________________________
-AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari()
-{
-  // constructor
-  fSegmentation=0;
-  fDigits=0;
-  fClusters=0;
-  fNclusters=0;
-  SetDx();
-  SetDz();
-  
-}
-
-//__________________________________________________________________________
-AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari(const
-AliITSClusterFinderSPDbari &source){
-  //     Copy Constructor 
-  if(&source == this) return;
-  this->fClusters = source.fClusters ;
-  this->fNclusters = source.fNclusters ;
-  this->fDz = source.fDz ;
-  this->fDx = source.fDx ;
-  return;
-}
-
-//_________________________________________________________________________
-AliITSClusterFinderSPDbari& 
-  AliITSClusterFinderSPDbari::operator=(const AliITSClusterFinderSPDbari &source) {
-  //    Assignment operator
-  if(&source == this) return *this;
-  this->fClusters = source.fClusters ;
-  this->fNclusters = source.fNclusters ;
-  this->fDz = source.fDz ;
-  this->fDx = source.fDx ;
-  return *this;
-}
-
-//_____________________________________________________________________________
-void AliITSClusterFinderSPDbari::FindRawClusters(Int_t module){
-   
-    // input of Cluster Finder  
-    Int_t digitcount=0;
-    Int_t numberd=100000;
-    Int_t   *digx       = new Int_t[numberd];
-    Int_t   *digz       = new Int_t[numberd];
-    Int_t   *digtr1     = new Int_t[numberd];
-    Int_t   *digtr2     = new Int_t[numberd];
-    Int_t   *digtr3     = new Int_t[numberd];
-    Int_t   *digtr4     = new Int_t[numberd];
-    
-    //  output of Cluster Finder    
-    Int_t numberc=10000;
-    Float_t *xcenterl   = new Float_t[numberc];
-    Float_t *zcenterl   = new Float_t[numberc];
-    Float_t *errxcenter = new Float_t[numberc];
-    Float_t *errzcenter = new Float_t[numberc];
-    Int_t   *tr1clus    = new Int_t[numberc];
-    Int_t   *tr2clus    = new Int_t[numberc];
-    Int_t   *tr3clus    = new Int_t[numberc];
-
-    Int_t nclus;
-
-    digitcount=0;
-    Int_t ndigits = fDigits->GetEntriesFast();  
-    if (!ndigits) return;
-
-
-    AliITSdigit *dig;
-    AliITSdigitSPD *dig1;
-    Int_t ndig;
-    for(ndig=0; ndig<ndigits; ndig++) {
-         dig= (AliITSdigit*)fDigits->UncheckedAt(ndig);
-     
-         digx[digitcount] = dig->fCoord2+1;  //starts at 1
-         digz[digitcount] = dig->fCoord1+1;  //starts at 1
-
-         dig1= (AliITSdigitSPD*)fDigits->UncheckedAt(ndig);
-
-         digtr1[digitcount] = dig1->fTracks[0];
-         digtr2[digitcount] = dig1->fTracks[1];
-         digtr3[digitcount] = dig1->fTracks[2];
-         digtr4[digitcount] = dig1->fSignal;
-
-         digitcount++;
-    }
-
-
-    ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,digtr4,
-              nclus,xcenterl,zcenterl,errxcenter,errzcenter,
-              tr1clus, tr2clus, tr3clus);
-    DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
-              tr1clus, tr2clus, tr3clus);
-
-
-    delete[] digx       ;
-    delete[] digz       ;
-    delete[] digtr1     ;
-    delete[] digtr2     ;
-    delete[] digtr3     ;
-    delete[] digtr4     ;
-    delete[] xcenterl   ;
-    delete[] zcenterl   ;
-    delete[] errxcenter ;
-    delete[] errzcenter ;
-    delete[] tr1clus    ;
-    delete[] tr2clus    ;
-    delete[] tr3clus    ;
-
-}
-//-----------------------------------------------------------------
-void AliITSClusterFinderSPDbari::ClusterFinder(Int_t ndigits,
-                    Int_t digx[],Int_t digz[],
-                    Int_t digtr1[],Int_t digtr2[],Int_t digtr3[],Int_t digtr4[],
-                    Int_t &nclus, Float_t xcenter[],Float_t zcenter[],
-                                   Float_t errxcenter[],Float_t errzcenter[],
-                    Int_t tr1clus[], Int_t tr2clus[], Int_t tr3clus[]) {
-//
-// Search for clusters of fired pixels (digits). Two digits are linked
-// inside a cluster if they are countiguous both in row or column
-// direction.  Diagonal digits are not linked.
-// xcenter, ycenter, zcenter are the coordinates of the center
-// of each found cluster, calculated from the averaging the corresponding
-// coordinate of the center of the linked digits. The coordinates are
-// given in the local reference sistem. 
-// errxcenter, errycenter, errzcenter are the errors associated to
-// the corresponding average.
-//
-//
-
-  Int_t if1, min, max, nd;
-  Int_t x1, z1, t1, t2, t3, t4;
-  Int_t ndx, ndz, ndxmin, ndxmax, ndzmin, ndzmax;
-  Float_t dx, dz; 
-  Int_t i,k,ipos=0;
-  Float_t xdum, zdum;      
-  Int_t kmax, sigmax;
-  Float_t deltax, deltaz;
-  Float_t ndig;
-
-  Int_t numberd=10000;
-  Int_t *ifpad = new Int_t[numberd];
-  Int_t *xpad  = new Int_t[numberd];
-  Int_t *zpad  = new Int_t[numberd];
-  Int_t *tr1pad  = new Int_t[numberd];
-  Int_t *tr2pad  = new Int_t[numberd];
-  Int_t *tr3pad  = new Int_t[numberd];
-  Int_t *tr4pad  = new Int_t[numberd];
-  Int_t *iclus   = new Int_t[numberd];
-
-  nclus=1;
-  for (i=0; i < ndigits ; i++){
-    ifpad[i] = -1;
-    iclus[i] = 0;
-  }
-
-  ifpad[0]=0;
-  for (i=0; i < ndigits-1 ; i++) {
-    if ( ifpad[i] == -1 ) { 
-               nclus++;
-               ipos++;
-       ifpad[i]=nclus-1;
-    }
-    for (Int_t j=i+1 ; j < ndigits ; j++)  {  
-      if (ifpad[j]== -1 ) {
-            dx = TMath::Abs(digx[i]-digx[j]);
-            dz = TMath::Abs(digz[i]-digz[j]);
-
-//          if ( ( dx+dz )==1 )  //clusters are not diagonal
-            if ( ( dx+dz )==1 || (dx==1 && dz==1) ) { //diagonal clusters allowed
-                   ipos++;
-                   ifpad[j]=ifpad[i];
-
-                   x1=digx[j];
-                   z1=digz[j];
-                   digx[j]=digx[ipos];
-                   digz[j]=digz[ipos];
-                   digx[ipos]=x1;
-                   digz[ipos]=z1;
-
-                   t1=digtr1[j];
-                   t2=digtr2[j];
-                   t3=digtr3[j];
-                   t4=digtr4[j];
-                   digtr1[j]=digtr1[ipos];
-                   digtr2[j]=digtr2[ipos];
-                   digtr3[j]=digtr3[ipos];
-                   digtr4[j]=digtr4[ipos];
-                   digtr1[ipos]=t1;
-                   digtr2[ipos]=t2;
-                   digtr3[ipos]=t3;
-                   digtr4[ipos]=t4;
-
-                   if1=ifpad[j];
-                   ifpad[j]=ifpad[ipos];
-                   ifpad[ipos]=if1;
-            }
-      }
-    }
-  }//end loop on digits   
-
-  if ( ifpad[ndigits-1] == -1 ) {
-         nclus++;
-         ifpad[ndigits-1]=nclus-1;
-  }
-
-  for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
-
-  min=0;
-  max=0;
-  // loop on found clusters 
-  for (i=0 ; i < nclus ; i++)  
-  {
-     min = max;
-     max += iclus[i];
-     deltax = fSegmentation->Dpx(0);
-     if (iclus[i]!=1) 
-     {
-        //cluster with more than one digit
-        nd=iclus[i];
-        ndig=(Float_t) nd;
-           Int_t count=0;
-        for (k=min;k<min+nd;k++)
-        {
-              xpad[count] = digx[k];      
-              zpad[count] = digz[k];
-
-              tr1pad[count] = digtr1[k];          
-              tr2pad[count] = digtr2[k];          
-              tr3pad[count] = digtr3[k];          
-              tr4pad[count] = digtr4[k];          
-
-              count++; 
-        }
-        ndxmin = xpad[TMath::LocMin(nd,xpad)];
-        ndxmax = xpad[TMath::LocMax(nd,xpad)];
-        ndzmin = zpad[TMath::LocMin(nd,zpad)];
-        ndzmax = zpad[TMath::LocMax(nd,zpad)];
-        ndx = ndxmax - ndxmin+1;
-        ndz = ndzmax - ndzmin+1;
-
-
-        // calculate x and z coordinates of the center of the cluster
-        fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
-
-           if (ndx == 1) {         
-                xcenter[i] = xdum;
-           }    
-           else{ 
-             xcenter[i] = 0.;
-                for (k=0;k<nd;k++) {
-                fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
-                   xcenter[i] += (xdum / nd);
-                }                     
-           }
-
-           if (ndz == 1) {
-                zcenter[i] = zdum;
-           }   
-           else {
-                zcenter[i] = 0.;
-                for (k=0;k<nd;k++) {         
-               fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
-                  zcenter[i] += (zdum / nd);
-                }
-           }
-
-        // error on points in x and z directions
-
-        if (ndx == 1) {
-             errxcenter[i] = deltax / TMath::Sqrt(12.);
-        }
-        else {
-             errxcenter[i] = 0.;                       
-             for (k=0;k<nd;k++){ 
-               fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
-               errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/(nd*(nd-1)); 
-             }   
-                errxcenter[i] = TMath::Sqrt(errxcenter[i]);
-        }
-       
-           if (ndz == 1) {
-            deltaz = fSegmentation->Dpz(digz[min]);                  
-               errzcenter[i] = deltaz / TMath::Sqrt(12.);
-        }
-           else {
-               errzcenter[i] = 0.;
-               for (k=0;k<nd;k++){ 
-               fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
-                  errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/(nd*(nd-1));
-               }
-               errzcenter[i] = TMath::Sqrt(errzcenter[i]);
-           }    
-        // take three track numbers for the cluster
-        // choose the track numbers of the digit with higher signal 
-        kmax = 0;
-        sigmax = 0;
-        for (k=0;k<nd;k++){
-          if(tr4pad[k] > sigmax){
-            sigmax = tr4pad[k];
-            kmax   = k;
-          }
-        }
-        if(sigmax != 0) {
-            tr1clus[i]= tr1pad[kmax];
-            tr2clus[i]= tr2pad[kmax];
-            tr3clus[i]= tr3pad[kmax];
-         }
-         else {
-            tr1clus[i]= -2;
-            tr2clus[i]= -2;
-            tr3clus[i]= -2;
-        }
-     }
-     else  {
-      
-        // cluster with single digit
-        ndig= 1.;
-           ndx = 1;
-        ndz = 1;
-        fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum,zdum);
-        xcenter[i] = xdum;
-        zcenter[i] = zdum;
-        tr1clus[i]=digtr1[min];
-        tr2clus[i]=digtr2[min];
-        tr3clus[i]=digtr3[min];
-           deltaz = fSegmentation->Dpz(digz[min]);
-           errxcenter[i] = deltax / TMath::Sqrt(12.);
-           errzcenter[i] = deltaz / TMath::Sqrt(12.);
-     }
-
-     // store the cluster information to the AliITSRawCLusterSPD object
-     AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
-
-     //put the cluster center in local reference frame of the detector
-     // and in microns
-     xcenter[i] = xcenter[i] - fSegmentation->Dx()/2.; 
-     zcenter[i] = zcenter[i] - fSegmentation->Dz()/2.;
-
-
-     AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i],xcenter[i],ndig,ndz,ndx,0.,0.,0.,0.,0.,0.,0.);
-     iTS->AddCluster(0,clust);
-     delete clust;
-  }//end loop on clusters   
-  delete[] ifpad;
-  delete[] xpad ;
-  delete[] zpad ;
-  delete[] iclus;
-  delete[] tr1pad;
-  delete[] tr2pad;
-  delete[] tr3pad;
-  delete[] tr4pad;
-}
-//______________________________________________________
-void AliITSClusterFinderSPDbari::DigitToPoint(Int_t nclus,
-                   Float_t *xcenter,   Float_t *zcenter,
-                  Float_t *errxcenter,Float_t *errzcenter, 
-                  Int_t *tr1clus, Int_t *tr2clus, Int_t *tr3clus){
- //
- // A point is associated to each cluster of SPD digits. The points
- // and their associated errors are stored in the file galiceSP.root.
- //
- //
-     Float_t l[3],xg,zg;
-     const Float_t kconv = 1.0e-4; // micron -> cm
-
-     // get rec points
-     AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
-
-     for (Int_t i=0; i<nclus; i++)
-     {
-        l[0] = kconv*xcenter[i];
-        l[1] = kconv*fSegmentation->Dy()/2.;
-        l[2] = kconv*zcenter[i];
-
-        xg = l[0]; 
-        zg = l[2]; 
-
-           Float_t sigma2x = (kconv*errxcenter[i]) * (kconv*errxcenter[i]);
-           Float_t sigma2z = (kconv*errzcenter[i]) * (kconv*errzcenter[i]);
-        AliITSRecPoint rnew;
-        rnew.SetX(xg);
-        rnew.SetZ(zg);
-        rnew.SetQ(1.);
-        rnew.SetdEdX(0.);
-        rnew.SetSigmaX2(sigma2x);
-        rnew.SetSigmaZ2(sigma2z);
-        rnew.fTracks[0]=tr1clus[i];
-        rnew.fTracks[1]=tr2clus[i];
-        rnew.fTracks[2]=tr3clus[i];
-        iTS->AddRecPoint(rnew); 
-     }
-}
diff --git a/ITS/AliITSClusterFinderSPDbari.h b/ITS/AliITSClusterFinderSPDbari.h
deleted file mode 100644 (file)
index cdf7c16..0000000
+++ /dev/null
@@ -1,64 +0,0 @@
-#ifndef ALIITSCLUSTERFINDERSPDBARI_H
-#define ALIITSCLUSTERFINDERSPDBARI_H
-
-////////////////////////////////////////////////
-//  ITS Cluster Finder Class                 //
-////////////////////////////////////////////////
-
-#include "AliITSClusterFinder.h"
-
-class AliITSMapA1;
-
-class AliITSClusterFinderSPDbari :
-  public AliITSClusterFinder
-
-{
-public:
-  AliITSClusterFinderSPDbari
-  (AliITSsegmentation *segmentation,
-   TClonesArray *digits, TClonesArray *recpoints);
-  AliITSClusterFinderSPDbari();
-  virtual ~AliITSClusterFinderSPDbari(){
-    // destructor
-  }
-  AliITSClusterFinderSPDbari(const AliITSClusterFinderSPDbari &source); // copy constructor
-  AliITSClusterFinderSPDbari& operator=(const AliITSClusterFinderSPDbari &source); // assignment operator
-  
-  
-  virtual void SetDx(Float_t dx=1.) {
-    // set dx
-    fDx=dx;
-  }
-  virtual void SetDz(Float_t dz=0.) {
-    // set dz
-    fDz=dz;
-  }
-
-  // Search for clusters
-  virtual void FindRawClusters(Int_t module); 
-  void  DigitToPoint(Int_t nclus, Float_t *xcenter, Float_t *zcenter,
-                  Float_t *errxcenter,Float_t *errzcenter,
-                 Int_t *tr1clus, Int_t *tr2clus, Int_t *tr3clus);
-  void  ClusterFinder(Int_t ndigits,Int_t digx[],Int_t digz[],
-                     Int_t digtr1[],Int_t digtr2[],Int_t digtr3[],Int_t digtr4[],
-                     Int_t &nclus,
-                     Float_t xcenter[],Float_t zcenter[],
-                     Float_t errxcenter[],Float_t errzcenter[],  
-                     Int_t tr1clus[],Int_t tr2clus[], Int_t tr3clus[]);  
-  
-  
-  
-private:
-  
-  TClonesArray       *fClusters;      // clusters
-  Int_t               fNclusters;     // num of clusters
-  Float_t             fDz;            // dz
-  Float_t             fDx;            // dx
-  
-  Int_t               fMinNCells;     // min num of cells in the cluster
-  
-  ClassDef(AliITSClusterFinderSPDbari,1)  // SPD clustering based
-                                          // on Nico Di Bari algorithm
-    };
-#endif
-
diff --git a/ITS/AliITSFindClustersBari.C b/ITS/AliITSFindClustersBari.C
deleted file mode 100644 (file)
index 7bbd301..0000000
+++ /dev/null
@@ -1,105 +0,0 @@
-Int_t AliITSFindClustersBari() {
-
-  printf("FindClusters\n");
-
-  TFile *in=TFile::Open("galice.root","UPDATE");
-  if (!in->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;}
-
-  in->ls();
-
-   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
-     cerr<<"gAlice have not been found on galice.root !\n";
-     return 2;
-   }
-
-
-   gAlice->GetEvent(0);
-
-   AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS"); 
-   if (!ITS) {
-     cerr<<"ITSFindClusters.C : AliITS object not found on file\n";
-     return 3;
-   }
-   Int_t ver = ITS->IsVersion(); 
-   cerr<<"ITS version "<<ver<<" has been found !\n";
-
-    ITS->MakeTreeC();
-// Set the models for cluster finding
-   AliITSgeom *geom = ITS->GetITSgeom();
-
-   // SPD
-   AliITSDetType *iDetType=ITS->DetType(0);
-   AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
-   TClonesArray *dig0  = ITS->DigitsAddress(0);
-   TClonesArray *recp0  = ITS->ClustersAddress(0);
-   AliITSClusterFinderSPDbari *rec0=new AliITSClusterFinderSPDbari(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);
-   AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
-   if (!seg1) seg1 = new AliITSsegmentationSDD(geom);
-   AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
-   if (!res1) res1=new AliITSresponseSDD();
-   Float_t baseline,noise;
-   res1->GetNoiseParam(noise,baseline);
-   Float_t noise_after_el = res1->GetNoiseAfterElectronics();
-   Float_t thres = baseline;
-   thres += (4.*noise_after_el);  // TB // (4.*noise_after_el);
-   printf("thres %f\n",thres);
-   res1->Print();
-   TClonesArray *dig1  = ITS->DigitsAddress(1);
-   TClonesArray *recp1  = ITS->ClustersAddress(1);
-   AliITSClusterFinderSDD *rec1=new AliITSClusterFinderSDD(seg1,res1,dig1,recp1);
-   rec1->SetCutAmplitude((int)thres);
-   ITS->SetReconstructionModel(1,rec1);
-
-
-   // SSD
-   AliITSDetType *iDetType=ITS->DetType(2);
-   AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
-   TClonesArray *dig2  = ITS->DigitsAddress(2);
-   AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
-   ITS->SetReconstructionModel(2,rec2);
-   // test
-   printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz());
-   printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
-
-
-
-   if(!gAlice->TreeR()) gAlice->MakeTree("R");
-   //make branch
-   ITS->MakeBranch("R");
-
-   TStopwatch timer;
-
-   switch (ver) {
-   case 5:
-      cerr<<"Looking for clusters...\n";
-      {
-       timer.Start();
-       ITS->DigitsToRecPoints(0,0,"All");
-      }
-      break;
-   default:
-      cerr<<"Invalid ITS version !\n";
-      return 5;
-   }
-
-   timer.Stop(); timer.Print();
-
-   delete rec0;
-   delete rec1;
-   delete rec2;
-
-
-   delete gAlice; gAlice=0;
-
-   in->Close();
-
-   return 0;
-}
diff --git a/ITS/AliITSHits2DigitsBari.C b/ITS/AliITSHits2DigitsBari.C
deleted file mode 100644 (file)
index 68cc01b..0000000
+++ /dev/null
@@ -1,114 +0,0 @@
-Int_t AliITSHits2DigitsBari()
-{
-
-  // Connect the Root Galice file containing Geometry, Kine and Hits
-
-  const char * inFile = "galice.root";  
-  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(inFile);
-  if (file) {file->Close(); delete file;}
-  printf("Hits2Digits\n");
-  file = new TFile(inFile,"UPDATE");
-  if (!file->IsOpen()) {
-    cerr<<"Can't open "<<inFile<<" !\n";
-    return 1;
-  }
-  file->ls();
-
-  // Get AliRun object from file or return if not on file
-  if (gAlice) delete gAlice;
-  gAlice = (AliRun*)file->Get("gAlice");
-  if (!gAlice) {
-    cerr<<"ITSHits2Digits.C : AliRun object not found on file\n";
-    return 2;
-  }
-
-  gAlice->GetEvent(0);
-  AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");      
-  if (!ITS) {
-    cerr<<"ITSHits2Digits.C : AliITS object not found on file\n";
-    return 3;
-  }
-
-// Set the simulation models for the three detector types
-  AliITSgeom *geom = ITS->GetITSgeom();
-
-  // SPD
-  AliITSDetType *iDetType=ITS->DetType(0);
-  AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
-  AliITSresponseSPDbari *res0 = (AliITSresponseSPDbari*)iDetType->GetResponseModel();
-  AliITSsimulationSPDbari *sim0=new AliITSsimulationSPDbari(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
-
-  // SDD
-  //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();
-  if (!res1) {
-    res1=new AliITSresponseSDD();
-    ITS->SetResponseModel(1,res1);
-  }
-   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
-  AliITSDetType *iDetType=ITS->DetType(2);
-  AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
-  AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
-  res2->SetSigmaSpread(3.,2.);
-  AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
-  ITS->SetSimulationModel(2,sim2);
-
-
-  cerr<<"Digitizing ITS...\n";
-  
-  TStopwatch timer;
-  timer.Start();
-  ITS->HitsToDigits(0,0,-1," ","All"," ");
-  timer.Stop(); timer.Print();
-
-  delete sim0;
-  delete sim1;
-  delete sim2;
-
-
-  delete gAlice;   gAlice=0;
-  file->Close(); 
-  delete file;
-  return 0;
-};
-
diff --git a/ITS/AliITSresponseSPDbari.cxx b/ITS/AliITSresponseSPDbari.cxx
deleted file mode 100644 (file)
index e49af35..0000000
+++ /dev/null
@@ -1,30 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-#include <TMath.h>
-
-#include "AliITSresponseSPDbari.h"
-//___________________________________________
-ClassImp(AliITSresponseSPDbari)        
-
-AliITSresponseSPDbari::AliITSresponseSPDbari()
-{
-  // constructor
-   SetThresholds();
-   SetNoiseParam();
-   SetDataType();
-}
-
diff --git a/ITS/AliITSresponseSPDbari.h b/ITS/AliITSresponseSPDbari.h
deleted file mode 100644 (file)
index f987d35..0000000
+++ /dev/null
@@ -1,70 +0,0 @@
-#ifndef ALIITSRESPONSESPDBARI_H
-#define ALIITSRESPONSESPDBARI_H
-
-#include "AliITSresponse.h"
-#include <TString.h>
-
-//----------------------------------------------
-//
-// ITS response class for SPD
-//
-class AliITSresponseSPDbari :
-  public AliITSresponse {
-public:
-  
-  AliITSresponseSPDbari();
-  virtual ~AliITSresponseSPDbari() { 
-    // destructror
-  }
-  //
-  // Configuration methods
-  //
-  
-  
-  virtual  void   SetDiffCoeff(Float_t p1=0) {
-    // 
-    fDiffCoeff=p1;
-  }
-  virtual  Float_t   DiffCoeff() {
-    // 
-    return fDiffCoeff;
-  }
-  virtual  void   SetThresholds(Float_t thresh=7.2e-6, Float_t sigma=1.0e-6) {
-    // Set Threshold and noise + threshold fluctuations parameter values
-    fThresh=thresh; fSigma=sigma;
-  }
-  virtual  void   Thresholds(Float_t &thresh, Float_t &sigma) {
-    // Get Threshold and noise + threshold fluctuations parameter values
-    thresh=fThresh; sigma=fSigma;
-  }
-  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() {
-    // Get data typer
-    return fDataType.Data();
-  } 
-  
-  ClassDef(AliITSresponseSPDbari,1) // SPD response
-    
-    protected:
-  
-  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
-};
-
-#endif
diff --git a/ITS/AliITSsimulationSPDbari.cxx b/ITS/AliITSsimulationSPDbari.cxx
deleted file mode 100644 (file)
index 1c3c3ce..0000000
+++ /dev/null
@@ -1,721 +0,0 @@
-#include <iostream.h>
-#include <TRandom.h>
-#include <TH1.h>
-#include <TMath.h>
-#include <TString.h>
-#include <TParticle.h>
-
-
-#include "AliRun.h"
-#include "AliITS.h"
-#include "AliITShit.h"
-#include "AliITSdigit.h"
-#include "AliITSmodule.h"
-#include "AliITSMapA2.h"
-#include "AliITSsimulationSPDbari.h"
-#include "AliITSsegmentation.h"
-#include "AliITSresponse.h"
-
-
-ClassImp(AliITSsimulationSPDbari)
-////////////////////////////////////////////////////////////////////////
-// Version: 0
-// Written by Rocco Caliandro
-// from a model developed with T. Virgili and R.A. Fini
-// June 15 2000
-//
-// AliITSsimulationSPD is the simulation of SPDs
-//
-//________________________________________________________________________
-
-AliITSsimulationSPDbari::AliITSsimulationSPDbari(){
-  // constructor
-  fResponse     = 0;
-  fSegmentation = 0;
-  fHis          = 0;
-  fThresh       = 0.;
-  fSigma        = 0.;
-  fCouplCol     = 0.;
-  fCouplRow     = 0.;
-}
-//_____________________________________________________________________________
-
-AliITSsimulationSPDbari::AliITSsimulationSPDbari(AliITSsegmentation *seg, AliITSresponse *resp) {
-  // constructor
-      fResponse = resp;
-      fSegmentation = seg;
-
-      fResponse->Thresholds(fThresh,fSigma);
-      fResponse->GetNoiseParam(fCouplCol,fCouplRow);
-      
-      fMapA2 = new AliITSMapA2(fSegmentation);
-   
-      // 
-      fNPixelsZ=fSegmentation->Npz();
-      fNPixelsX=fSegmentation->Npx();
-      fHis=0;
-}
-
-//_____________________________________________________________________________
-
-AliITSsimulationSPDbari::~AliITSsimulationSPDbari() { 
-  // destructor
-
-  delete fMapA2;
-
-  if (fHis) {
-     fHis->Delete(); 
-     delete fHis;     
-  }                
-}
-
-//__________________________________________________________________________
-AliITSsimulationSPDbari::AliITSsimulationSPDbari(const AliITSsimulationSPDbari &source){
-  //     Copy Constructor 
-  if(&source == this) return;
-  this->fMapA2 = source.fMapA2;
-  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;
-  return;
-}
-
-//_________________________________________________________________________
-AliITSsimulationSPDbari& 
-  AliITSsimulationSPDbari::operator=(const AliITSsimulationSPDbari &source) {
-  //    Assignment operator
-  if(&source == this) return *this;
-  this->fMapA2 = source.fMapA2;
-  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;
-  return *this;
-  }
-//_____________________________________________________________________________
-
-void AliITSsimulationSPDbari::DigitiseModule(AliITSmodule *mod, Int_t module,
-                                             Int_t dummy) {
-  // digitize module
-
-
-  TObjArray *fHits = mod->GetHits();
-  Int_t nhits = fHits->GetEntriesFast();
-  if (!nhits) return;
-
-
-  //printf("Module %d (%d hits) \n",module+1,nhits);
-
-
-  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];
-
-  // 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);
-
-
-
-  // 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 AliITSsimulationSPDbari::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 AliITSsimulationSPDbari::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)
-//
-
-
-   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
-
-   TObjArray *fHits = mod->GetHits();
-   AliITShit *hit = (AliITShit*) fHits->At(hitpos);
-   layer = hit->GetLayer();
-   etot=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
-    */
-    
-   
-
-        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;                 
-        }
-}
-
-//_________________________________________________________________________
-
-void AliITSsimulationSPDbari::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;
-
-
-   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);
-
-   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    
-
-   Float_t xpos, zpos;
-
-   fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos); 
-
-   Float_t xsize = fSegmentation->Dpx(0);
-   Float_t zsize = fSegmentation->Dpz(r1-1);
-
-   if (dirx == 1) refr = xpos+xsize/2.;
-             else refr = xpos-xsize/2.;
-
-   if (dirz == 1) refn = zpos+zsize/2.;
-             else refn = zpos-zsize/2.;
-
-   
-   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;
-
-      }
-      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);
-
-}
-//___________________________________________________________________________
-void AliITSsimulationSPDbari::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);          
-   
-   }
-
-}
-//___________________________________________________________________________
-void AliITSsimulationSPDbari::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;  
-                 signal = signal*1.0e9;  //signal in eV
-                 digits[2] =  (Int_t) signal;  // the signal is stored in eV
-                 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
-    }
-    
-}
-//_____________________________________________________________________________
-
-void AliITSsimulationSPDbari::GetList(Int_t label,Int_t idhit, Float_t **pList,
-                                      Int_t row, Int_t col) {
-  // loop over nonzero digits
-
-  Int_t ix = col;
-  Int_t iz = row;
-  Int_t globalIndex;
-  Float_t signal;
-  Float_t highest,middle,lowest;
-
-          
-  signal=fMapA2->GetSignal(iz,ix);
-
-
-  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];
-
-
-     // 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.;
-
-     *pList[globalIndex] = (float)label;
-     *(pList[globalIndex]+3) = signal;
-     *(pList[globalIndex]+6) = (float)idhit;
-  }
-  else{
-
-
-         // check the signal magnitude
-      highest = *(pList[globalIndex]+3);
-      middle  = *(pList[globalIndex]+4);
-      lowest  = *(pList[globalIndex]+5);
-
-
-      signal -= (highest+middle+lowest);
-
-
-         //
-         //  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 AliITSsimulationSPDbari::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 AliITSsimulationSPDbari::CreateHistograms() {
-  // CreateHistograms
-
-      Int_t i;
-      fHis=new TObjArray(fNPixelsZ);
-      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 AliITSsimulationSPDbari::ResetHistograms() {
-    //
-    // Reset histograms for this detector
-    //
-    Int_t i;
-    for(i=0;i<fNPixelsZ;i++ ) {
-       if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
-    }
-
-}
diff --git a/ITS/AliITSsimulationSPDbari.h b/ITS/AliITSsimulationSPDbari.h
deleted file mode 100644 (file)
index d36873a..0000000
+++ /dev/null
@@ -1,64 +0,0 @@
-#ifndef ALIITSSIMULATIONSPDBARI_H
-#define ALIITSSIMULATIONSPDBARI_H
-
-#include "AliITSsimulation.h"
-
-class AliITSMapA2;
-class AliITSsegmentation;
-class AliITSresponse;
-class AliITSmodule;
-
-//-------------------------------------------------------------------
-
-class AliITSsimulationSPDbari : public AliITSsimulation {
-
-public:
-        
-  AliITSsimulationSPDbari();
-  AliITSsimulationSPDbari(AliITSsegmentation *seg, AliITSresponse *res);
-  ~AliITSsimulationSPDbari();
-  AliITSsimulationSPDbari(const AliITSsimulationSPDbari &source); // copy constructor
-  AliITSsimulationSPDbari& operator=(const AliITSsimulationSPDbari &source); // ass. operator
-
-  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();
-  TObjArray*  GetHistArray() {
-    // get hist array
-    return fHis;
-  }
-
-private:
-
-  AliITSMapA2  *fMapA2;        // MapA2
-  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
-
-  TObjArray *fHis;             // just in case for histogramming
-    
-  ClassDef(AliITSsimulationSPDbari,1)  // Simulation of SPD clusters
-
-};
-
-#endif 
-
diff --git a/ITS/AliITStestBari.C b/ITS/AliITStestBari.C
deleted file mode 100644 (file)
index 07cd575..0000000
+++ /dev/null
@@ -1,35 +0,0 @@
-Int_t AliITStestBari() {
-   Int_t rc=0;
-
-//Test ITS simulation
-   gROOT->LoadMacro("$(ALICE_ROOT)/macros/grun.C");
-   grun();
-
-   Int_t ver=gAlice->GetDetector("ITS")->IsVersion();
-   delete gAlice; gAlice=0;
-
-   if (ver!=5) {
-      cerr<<"Invalid ITS version: "<<ver<<" ! (must be 5 for the moment)\n";
-      return 12345;
-   }
-
-   if (ver==5) {
-     gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2DigitsBari.C");
-     if (rc=AliITSHits2DigitsBari()) return rc;
-
-   }
-
-   printf("start reconstruction\n");
-
-//Test ITS reconstruction
-   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersBari.C");
-
-   delete gAlice; gAlice=0;
-   
-   if (rc=AliITSFindClustersBari()) return rc;
-
-   //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSgraphycs.C");
-   //if (rc=AliITSgraphycs()) return rc;
-
-   return rc;
-}
diff --git a/ITS/ITSDigitsToClustersBari.C b/ITS/ITSDigitsToClustersBari.C
deleted file mode 100644 (file)
index 9713586..0000000
+++ /dev/null
@@ -1,150 +0,0 @@
-#include "iostream.h"
-
-void ITSDigitsToClustersBari (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");
-   printf("file %p\n",file);
-   if (file) file->Close(); 
-   file = new TFile("galice.root","UPDATE");
-   file->ls();
-
-   printf ("I'm after Map \n");
-
-// 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");
-   }
-   printf ("I'm after gAlice \n");
-   
-   AliITS *ITS  = (AliITS*) gAlice->GetModule("ITS");
-   if (!ITS) return;
-
-   AliITSgeom *geom = ITS->GetITSgeom();
-
-
-   // NOTE: if you foresee to have (in segmentation or response) parameter
-   // values other than the default ones, and these values are used not only in
-   // 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 (ugly and impractical, no? )
-
-
-
-   // Set the models for cluster finding
-
-   // 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);
-   AliITSClusterFinderSPDbari *rec0=new AliITSClusterFinderSPDbari(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);
-   AliITSgeom *geom = ITS->GetITSgeom();
-
-   AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
-   if (!seg1) seg1 = new AliITSsegmentationSDD(geom);
-   AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
-   if (!res1) res1=new AliITSresponseSDD();
-
-   Float_t baseline,noise;
-   res1->GetNoiseParam(noise,baseline);
-   Float_t noise_after_el = res1->GetNoiseAfterElectronics();
-   Float_t thres = baseline;
-   thres += (4.*noise_after_el);  // TB // (4.*noise_after_el);
-   printf("thres %f\n",thres);
-   res1->Print();
-
-   TClonesArray *dig1  = ITS->DigitsAddress(1);
-   TClonesArray *recp1  = ITS->ClustersAddress(1);
-   AliITSClusterFinderSDD *rec1=new AliITSClusterFinderSDD(seg1,res1,dig1,recp1);
-   rec1->SetCutAmplitude((int)thres);
-   ITS->SetReconstructionModel(1,rec1);
-   rec1->Print();
-
-   // SSD
-
-   AliITSDetType *iDetType=ITS->DetType(2);
-   AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
-   TClonesArray *dig2  = ITS->DigitsAddress(2);
-   AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
-   ITS->SetReconstructionModel(2,rec2);
-   // test
-   //printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz());
-   //printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
-
-
-
-//
-// Event Loop
-//
-
-   cout << "Looking for clusters...\n";
-   
-   TStopwatch timer;
-
-   if(!gAlice->TreeR()) gAlice->MakeTree("R");
-   //make branch
-   ITS->MakeBranch("R");
-
-   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");
-       }     
-       cout << "nev         " <<nev<<endl;
-       cout << "nparticles  " <<nparticles<<endl;
-       if (nev < evNumber1) continue;
-       if (nparticles <= 0) return;
-
-       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 
-
-   delete rec0;
-   delete rec1;
-   delete rec2;
-
-   file->Close();
-}
diff --git a/ITS/ITSHitsToDigitsBari.C b/ITS/ITSHitsToDigitsBari.C
deleted file mode 100644 (file)
index ee01aed..0000000
+++ /dev/null
@@ -1,170 +0,0 @@
-#include "iostream.h"
-
-void ITSHitsToDigitsBari (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal  =25, Int_t size=-1) 
-{
-/////////////////////////////////////////////////////////////////////////
-//   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");
-   printf("file %p\n",file);
-   if (file) file->Close(); 
-   if (!file) file = new TFile("galice.root","UPDATE");
-   file->ls();
-
-   printf ("I'm after Map \n");
-
-// 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");
-   }
-   printf ("I'm after gAlice \n");
-   
-   AliITS *ITS  = (AliITS*) gAlice->GetModule("ITS");
-   if (!ITS) return;
-
-
-   // Set the simulation models
-
-   AliITSgeom *geom = ITS->GetITSgeom();
-
-   // SDD
-   // SDD compression param: 2 fDecrease, 2fTmin, 2fTmax or disable, 2 fTolerance
-                
-   AliITSDetType *iDetType=ITS->DetType(1);
-   AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
-   if (!res1) {
-         res1=new AliITSresponseSDD();
-         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;
-   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->SetElectronics(2);  // 1 = Pascal, 2 = OLA
-
-   res1->Print();
-
-   //cout << "SDD segmentation" << endl;
-
-   AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
-   if (!seg1) {
-       seg1 = new AliITSsegmentationSDD(geom,res1);
-       ITS->SetSegmentationModel(1,seg1);
-   }
-   seg1->Print();
-
-   //cout << "SDD segmentation" << endl;
-   AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1);
-   ITS->SetSimulationModel(1,sim1);
-   sim1->Print();
-
-
-   // SPD
-
-   AliITSDetType *iDetType=ITS->DetType(0);
-   AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
-   //AliITSresponseSPDbari *res0 = (AliITSresponseSPDbari*)iDetType->GetResponseModel();
-   AliITSresponseSPDbari *res0= new AliITSresponseSPDbari();
-   ITS->SetResponseModel(0,res0);
-
-// to change the  parameters
-   //res0->SetThresholds(7.2e-6, 1.e-6);
-   //res0->SetNoiseParam(0., 0.);
-   //res0->SetNoiseParam(0.04, 0.08);
-
-// to monitor the  parameters
-   Float_t thresh, sigma;
-   res0->Thresholds(thresh, sigma);
-   printf("SPDbari: threshold %e sigma  %e\n",thresh, sigma);
-   Float_t col, row;
-   res0->GetNoiseParam(col, row);
-   printf("SPDbari: Coupling by column %e Coupling by row %e\n",col, row);
-
-   AliITSsimulationSPDbari *sim0=new AliITSsimulationSPDbari(seg0,res0);
-   ITS->SetSimulationModel(0,sim0);
-
-   // SSD
-
-   AliITSDetType *iDetType=ITS->DetType(2);
-   AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
-   AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
-   res2->SetSigmaSpread(3.,2.);
-   AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
-   ITS->SetSimulationModel(2,sim2);
-
-
-//
-// Event Loop
-//
-
-
-   // create the TreeD 
-
-   Int_t nparticles=gAlice->GetEvent(0);
-   printf("Create TreeD \n");
-   if(!gAlice->TreeD()) gAlice->MakeTree("D");
-   //make branch
-   ITS->MakeBranch("D");
-
-   Int_t nbgr_ev=0;
-       
-       
-       cout<<"Digitizing ITS...\n";
-   TStopwatch timer;
-       
-   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");
-       }
-       cout << "nparticles  " <<nparticles<<endl;
-       if (nev < evNumber1) continue;
-       if (nparticles <= 0) return;
-
-       Int_t nbgr_ev=0;
-       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 sim1;
-   delete sim2;
-
-   file->Close();
-}
diff --git a/ITS/SPDclusterTestBari.C b/ITS/SPDclusterTestBari.C
deleted file mode 100644 (file)
index efaf903..0000000
+++ /dev/null
@@ -1,863 +0,0 @@
-void SPDclusterTestBari(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_bari.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;
- }
-
-   
- // 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]);
-    }
-   }
-
-    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();
-//
-   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();
-//
-   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();
-}
-