Updated version of the Bari code to work with the HEAD. A new test macros has also...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 23 May 2001 14:14:00 +0000 (14:14 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 23 May 2001 14:14:00 +0000 (14:14 +0000)
ITS/AliITSClusterFinderSPDbari.cxx
ITS/AliITSClusterFinderSPDbari.h
ITS/AliITSsimulationSPDbari.cxx
ITS/ITSReadPlotData.C
ITS/ITSsddanalysis.C
ITS/ITSsddtest.C
ITS/README
ITS/SPDclusterTestBari.C [new file with mode: 0644]
ITS/SSDrecpointTest.C

index efbb5d7..fc1209f 100644 (file)
@@ -80,7 +80,7 @@ AliITSClusterFinderSPDbari&
 }
 
 //_____________________________________________________________________________
-void AliITSClusterFinderSPDbari::FindRawClusters(){
+void AliITSClusterFinderSPDbari::FindRawClusters(Int_t module){
    
     // input of Cluster Finder  
     Int_t digitcount=0;
index 9992943..cdf7c16 100644 (file)
@@ -35,7 +35,7 @@ public:
   }
 
   // Search for clusters
-  virtual void FindRawClusters(); 
+  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);
index e102b77..8bc349c 100644 (file)
@@ -185,52 +185,52 @@ void AliITSsimulationSPDbari::HitToDigit(AliITSmodule *mod, Int_t hitpos, Int_t
    ntrack=hit->GetTrack();
    idhit=mod->GetHitHitIndex(hitpos);     
 
-   /* //debug
-     printf("layer,etot,ntrack,status %d %f %d %d\n",layer,etot,ntrack,hit->GetTrackStatus()); //debug
+    
+    /*
+    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 83750 effective sensitive area
-          // introduced in geometry (Dz=83600)
-          if (z1l>fSegmentation->Dz()) return;
-          if (z2l>fSegmentation->Dz()) return;
-          // to preserve for hit having x>12800
-          if (x1l>fSegmentation->Dx()) return;
-          if (x2l>fSegmentation->Dx()) return;
+          // 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;
 
-          // to account for tracks at the edge of the sensitive area
-          // of SPDs
-          if (x1l<0) return;
-          if (x2l<0) return;
-          if (z1l<0) return;
-          if (z2l<0) return;
-             
 
              if ((r1==r2) && (c1==c2)) 
              {
@@ -319,7 +319,7 @@ void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
    fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos); 
 
    Float_t xsize = fSegmentation->Dpx(0);
-   Float_t zsize = fSegmentation->Dpz(r1);
+   Float_t zsize = fSegmentation->Dpz(r1-1);
 
    if (dirx == 1) refr = xpos+xsize/2.;
              else refr = xpos-xsize/2.;
@@ -367,7 +367,10 @@ void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
             }     
 
          // shift to the pixel in the next cell in row direction
-         Float_t zsizeNext = fSegmentation->Dpz(rb);
+         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;
 
       }
@@ -387,7 +390,10 @@ void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
             }
 
          // shift to the pixel in the next cell in column direction
-         Float_t xsizeNext = fSegmentation->Dpx(cb);
+         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;
         
       }
@@ -397,8 +403,6 @@ void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
       epar = etot*(epar/dtot);
 
       //store row, column and energy lost in the crossed pixel
-
-
       frowpixel[npixel] = r1;
       fcolpixel[npixel] = c1;
       fenepixel[npixel] = epar;
@@ -564,6 +568,7 @@ void AliITSsimulationSPDbari::GetList(Int_t label,Int_t idhit, Float_t **pList,
           
   signal=fMapA2->GetSignal(iz,ix);
 
+
   globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 1
 
 
index c26ede8..9a382db 100644 (file)
@@ -179,6 +179,8 @@ Int_t ITSReadPlotData(char *filename = "galice.root", Int_t evNum = 0) {
                        cout << "No hits in module!" << endl;
                        continue;
                }
+               else
+                       cout << "Hits scanned..." << endl;
                for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]);
                
                // Reads recpoints...
@@ -187,6 +189,8 @@ Int_t ITSReadPlotData(char *filename = "galice.root", Int_t evNum = 0) {
                        cout << "No recpoints in module!" << endl;
                        continue;
                }
+               else
+                       cout << "Recpoints scanned..." << endl;
                for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]);
                
                // Reads digits...
@@ -195,6 +199,8 @@ Int_t ITSReadPlotData(char *filename = "galice.root", Int_t evNum = 0) {
                        cout << "No digits in module!" << endl;
                        //continue;
                }
+               else
+                       cout << "Digits scanned..." << endl;
                for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]);
 
                // Draws read data...
@@ -234,6 +240,22 @@ Int_t ITSReadPlotData(char *filename = "galice.root", Int_t evNum = 0) {
                legend->Draw();
                
                viewer->Update();
+               
+               
+               Text_t fname[250],ans;
+               cout << "Do you want to save the current canvas on a file (y/n) ? ";
+               cin >> ans;
+               if(ans == 'y' || ans == 'Y') {
+                  cout << "Enter filename: ";
+                  cin >> fname;
+                  TString *control = new TString(fname);
+                  Bool_t ok=control->Contains(".C") || control->Contains(".root") || control->Contains(".ps") || control->Contains(".eps") || control->Contains(".gif");
+                  if(!ok){ 
+                     cout << "File extension is not recognized. The canvas will be saved as Postscript file";
+                     strcat(fname,".ps");
+                  }  
+                  viewer->SaveAs(fname);
+               }
        }
        
        cout << "Done. Goodbye" << endl;
@@ -332,6 +354,7 @@ Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*
        if (!digits_num)
                return 0;
        else {
+               cout << "Digits to scan: " << digits_num << endl;
                if (X) delete [] X;                     
                if (Z) delete [] Z;
                X = new Float_t[digits_num];            
@@ -346,13 +369,18 @@ Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*
                }
        }
   for (Int_t j = 0; j < digits_num; j++) {
-       cout << j << endl;
                digit = (AliITSdigit*)digits_array->UncheckedAt(j);
                Int_t iz=digit->fCoord1;  // cell number z
                Int_t ix=digit->fCoord2;  // cell number x
     // Get local coordinates of the element (microns)
-               if(dtype < 2)
-       seg->DetToLocal(ix, iz, X[j], Z[j]);
+               // ******************************* PARTE CORRETTA ***************************************
+               if(dtype < 2) {
+                       Float_t xx, zz; // aggiunta
+       seg->DetToLocal(ix, iz, xx, zz);
+                       X[j] = xx; // aggiunta
+                       Z[j] = zz; // aggiunta
+               }
+               // ******************************* FINE PARTE CORRETTA ***************************************
     else {
                        // SSD: if iz==0 ---> N side; if iz==1 P side
       if (ssdone[j] == 0) {
@@ -375,8 +403,8 @@ Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*
                                        }
                                }
         if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]);
-       X[j] /= 10000.0;  // convert microns to cm
-       Z[j] /= 10000.0;  // convert microns to cm
+                               X[j] /= 10000.0;  // convert microns to cm
+                               Z[j] /= 10000.0;  // convert microns to cm
                        }
                }
        }
index ed94d88..5db71d9 100644 (file)
@@ -225,19 +225,12 @@ void ITSsddanalysis (Int_t evNumber1=0,Int_t evNumber2=0)
     ITS->GetTreeC(nev);
     TTree *TC=ITS->TreeC();
     Int_t nentclu=TC->GetEntries();
-    printf("Found %d entries in the Hit tree (must be one per module per event!)\n",nentclu);
+    printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);
     
     TTree *TR = gAlice->TreeR();
     Int_t nentrec=TR->GetEntries();
     printf("Found %d entries in the Rec tree (must be one per module per event!)\n",nentrec);
     
-    // Get Pointers to Clusters & Recpoints TClonesArrays
-    
-    Int_t iDet = 1;  // 1 = SDD
-
-    TClonesArray *ITSclu  = ITS->ClustersAddress(iDet);
-    TClonesArray *ITSrec  = ITS->RecPoints(); 
-    
     // check recpoints
     
     Int_t nbytes = 0;
@@ -252,12 +245,25 @@ void ITSsddanalysis (Int_t evNumber1=0,Int_t evNumber2=0)
     
     TObjArray *fITSmodules = ITS->GetModules();
     
-    Int_t first = aliitsgeo->GetStartDet(iDet);
-    Int_t last = aliitsgeo->GetLastDet(iDet);
+    Int_t iDet = 1;  // 1 = SDD
+
+    Int_t first = aliitsgeo->GetStartDet(0);
+    Int_t last = aliitsgeo->GetLastDet(2);
+    Int_t first_ok = aliitsgeo->GetStartDet(iDet);
+    Int_t last_ok = aliitsgeo->GetLastDet(iDet);
     printf("det type %d first, last %d %d \n",iDet,first,last);
     
-    for (Int_t mod=0; mod<last-first+1; mod++) {
-      cout << "Module: " << mod+1 << endl;
+    TClonesArray *ITSclu = 0;
+    TClonesArray *ITSrec = 0;
+    for (Int_t mod=first_ok; mod<=last_ok; mod++) {
+      cout << "Module: " << mod << endl;
+      //      if(mod < first_ok) continue;
+      //      if(mod > last_ok) continue;
+      // Get Pointers to Clusters & Recpoints TClonesArrays
+    
+      ITSclu  = ITS->ClustersAddress(iDet);
+      ITSrec  = ITS->RecPoints(); 
+    
       TTree *TR = gAlice->TreeR();
       Int_t nentrec=TR->GetEntries();
       TClonesArray *ITSrec  = ITS->RecPoints(); 
@@ -269,12 +275,12 @@ void ITSsddanalysis (Int_t evNumber1=0,Int_t evNumber2=0)
       
       Int_t nrecp = ITSrec->GetEntries();
       totpoints += nrecp;
-      //if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
+      if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
       if (!nrecp) continue;
       
       Int_t nrecc = ITSclu->GetEntries();
       totclust += nrecc;
-      //if (nrecc) printf("Found %d clusters for module %d\n",nrecc,mod);
+      if (nrecc) printf("Found %d clusters for module %d\n",nrecc,mod);
       
       Int_t nrecp = ITSrec->GetEntries();
       Int_t startSDD = aliitsgeo->GetStartSDD();
index 8e32619..3aad476 100644 (file)
@@ -2,16 +2,6 @@
 
 gROOT->Reset();
 
-// Dynamically link some shared libs
-
-  if (gClassTable->GetID("AliRun") < 0) {
-    gROOT->LoadMacro("loadlibs.C");
-    loadlibs();
-  } else {
-      delete gAlice;
-      gAlice=0;
-  }
-
 TFile f("SDD_histos_test.root");
 
 Int_t nbins = 18;
index 63218a3..d830564 100644 (file)
@@ -3,19 +3,20 @@
 
 Dear ITS (and ALICE) user,
 
-this is a short help on how to run the ITS simulation/reconstruction code 
+     This is a short help on how to run the ITS simulation/reconstruction code 
 within the AliRoot framework. It is NOT intended as a comprehensive user's 
 guide and eventually it will be updated in the ITS Manual which is on its way 
-to be written. What follows requires that the you already know how to download, 
+to be written. What follows requires that the you already know how to download,
 install and compile the AliRoot package. This file explains how to set the 
 proper ITS geometry and how to tun the test macros contained under the 
 directory ITS in order to compare your own installation with the standard one 
-looking at some distributions/histograms. 
-Any difference between what is described here and what you really get when you 
-run the code on your computer must be reported to Roberto Barbera at 
-roberto.barbera@ct.infn.it. Please note that all the tests described here have 
-been done on a PC running Linux RedHat 6.1, gcc 2.95.2, and Root 3.00/06. 
-If you have different hardware/software configuration, please add it to all bug
+looking at some distributions/histograms.
+        Any difference between what is described here and what you really get
+when you  run the code on your computer must be reported either to Roberto
+Barbera at roberto.barbera@ct.infn.it or to Bjorn Nilsen at
+nilsen@mps.ohio-state.edu. Please note that all the tests described here have
+been done on a PC running Linux RedHat 6.1, gcc 2.95.2, and Root 3.00/06. If
+you have different hardware/software configuration, please add it to all bug
 reports.
 
 
@@ -23,28 +24,27 @@ reports.
 ---------------------------------------
 
 In order to set one of the many ITS geometries available, you have to modify 
-the ITS part in the file Config.C under the directory macros (you have to modify
-the file ConfigPPR.C if you want to run full 'PPR' events). The default ITS 
-part of Config.C (or ConfigPPR.C) is reported here:
+the ITS part in the file Config.C under the directory macros (you have to
+modify the file ConfigPPR.C if you want to run full 'PPR' events). The default
+ITS part of Config.C (or ConfigPPR.C) is reported here:
 
   if(iITS) {
 
 //=================== ITS parameters ============================
     //
-    // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
-    // almost all other detectors. This involves the fact that the ITS geometry
-    // still has several options to be followed in parallel in order to determine
-    // the best set-up which minimizes the induced background. All the geometries
-    // available to date are described in the following. Read carefully the comments
-    // and use the default version (the only one uncommented) unless you are making
-    // comparisons and you know what you are doing. In this case just uncomment the
-    // ITS geometry you want to use and run Aliroot.
+    // As the innermost detector in ALICE, the Inner Tracking System "impacts"
+    // on almost all other detectors. This involves the fact that the ITS
+    // geometry  still has several options to be followed in parallel in order
+    // to determine the best set-up which minimizes the induced background. All
+    // the geometries available to date are described in the following. Read
+    // carefully the comments and use the default version (the only one
+    // uncommented) unless you are making comparisons and you know what you are
+    // doing. In this case just uncomment the ITS geometry you want to use and
+    // run Aliroot.
     //
     // Detailed geometries:         
     //
     //
-    //
-    //
     //AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
     //
     AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
@@ -57,8 +57,8 @@ part of Config.C (or ConfigPPR.C) is reported here:
     //ITS->SetThicknessDet2(300.);   // detector thickness on layer 2 must be in the range [100,300]
     //ITS->SetThicknessChip1(300.);  // chip thickness on layer 1 must be in the range [150,300]
     //ITS->SetThicknessChip2(300.);  // chip thickness on layer 2 must be in the range [150,300]
-    //ITS->SetRails(1);                   // 1 --> rails in ; 0 --> rails out
-    //ITS->SetCoolingFluid(1);    // 1 --> water ; 0 --> freon
+    //ITS->SetRails(1);            // 1 --> rails in ; 0 --> rails out
+    //ITS->SetCoolingFluid(1);     // 1 --> water ; 0 --> freon
     //
     //AliITSvPPRsymm *ITS  = new AliITSvPPRsymm("ITS","New ITS PPR detailed version with symmetric services");
     //ITS->SetMinorVersion(2);                                      
@@ -76,7 +76,6 @@ part of Config.C (or ConfigPPR.C) is reported here:
     // for reconstruction !):
     //                                                     
     //
-    //
     //AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS coarse version with asymmetric services");
     //ITS->SetRails(1);                // 1 --> rails in ; 0 --> rails out
     //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
@@ -90,8 +89,8 @@ part of Config.C (or ConfigPPR.C) is reported here:
     // Geant3 <-> EUCLID conversion
     // ============================
     //
-    // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
-    // media to two ASCII files (called by default ITSgeometry.euc and
+    // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry
+    // and  media to two ASCII files (called by default ITSgeometry.euc and
     // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
     // The default (=0) means that you dont want to use this facility.
     //
@@ -100,62 +99,109 @@ part of Config.C (or ConfigPPR.C) is reported here:
   
 As you can see looking at the uncommented line, the present default is
 
-    AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
+    AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with
+asymmetric services"); 
     
-which is the TDR detailed geometry with asymmetric services. 
-If you want to run the TDR detailed version with symmetric services, the only 
-uncommented line must be:
+which is the TDR detailed geometry with asymmetric services. If you want to run
+the TDR detailed version with symmetric services, You have to uncomment this
+line and comment out the line above.
 
-    AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
+    AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with
+symmetric services");
 
-If you want to run the new PPR coarse geometry with asymmetric services, the 
-only uncommented lines must be:
+If you want to run the new PPR coarse geometry with asymmetric services, You
+have to uncomment this line and comment out the line above.
 
-    AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS coarse version with asymmetric services");
+    AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS coarse
+version with asymmetric services");
     ITS->SetRails(1);                // 1 --> rails in ; 0 --> rails out
     ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
 
-If you want to run the new PPR coarse geometry with symmetric services, the 
-only uncommented lines must be:
-
-    AliITS *ITS  = new AliITSvPPRcoarsesymm("ITS","New ITS coarse version with symmetric services");
+There are two functions shown above. The first set the ITS rails in (=1) or out
+(=0) (the default is in). If the switch for rails is different from 0 and 1 a
+warning message is printed out and the default (rails in) is used. The second
+changes the material of the supports to the ITS services': copper (=0, which is
+the default), aluminum (=1), and  carbon (=2). If the switch of the support
+material is different from 0, 1 and 2, a warning message is printed out and the
+default (copper) is used. Note the the thickness of the supports is always the
+same so the situation is completely unrealistic. This is because there has been
+no attempt to change the geometry of this material as would be needed for the
+support to realy support the serveces. The possibility to  change this material
+has been explicitly requested by the PMD group to allow them to study the
+effect of this material on their detector.
+
+If you want to run the new PPR coarse geometry with symmetric services, you
+have to uncomment these lines and comment out the above lines.
+
+    AliITS *ITS  = new AliITSvPPRcoarsesymm("ITS","New ITS coarse version with
+symmetric services");
     ITS->SetRails(1);                // 1 --> rails in ; 0 --> rails out
     ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
 
-If you want to run the new PPR detailed geometry with asymmetric services, the 
-only uncommented lines must be:
+The two additional functions are just the same of those described above.
+
+    If you want to run the new PPR detailed geometry with asymmetric services,
+you have to uncomment these lines and comment out the above lines.
 
-    AliITSvPPRasymm *ITS  = new AliITSvPPRasymm("ITS","New ITS PPR detailed version with asymmetric services");
+    AliITSvPPRasymm *ITS  = new AliITSvPPRasymm("ITS","New ITS PPR detailed version
+with asymmetric services");
     ITS->SetMinorVersion(2);
     ITS->SetReadDet(kFALSE);
     ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det");
-    ITS->SetThicknessDet1(300.);   // detector thickness on layer 1 must be in the range [100,300]
-    ITS->SetThicknessDet2(300.);   // detector thickness on layer 2 must be in the range [100,300]
-    ITS->SetThicknessChip1(300.);  // chip thickness on layer 1 must be in the range [150,300]
-    ITS->SetThicknessChip2(300.);  // chip thickness on layer 2 must be in the range [150,300]
-    ITS->SetRails(1);             // 1 --> rails in ; 0 --> rails out
-    ITS->SetCoolingFluid(1);      // 1 --> water ; 0 --> freon
-
-If you want to run the new PPR detailed geometry with symmetric services, the 
-only uncommented lines must be:
+    ITS->SetThicknessDet1(300.);   // detector thickness on layer 1 must be in the
+range [100,300]
+    ITS->SetThicknessDet2(300.);   // detector thickness on layer 2 must be in the
+range [100,300]
+    ITS->SetThicknessChip1(300.);  // chip thickness on layer 1 must be in the range
+[150,300]
+    ITS->SetThicknessChip2(300.);  // chip thickness on layer 2 must be in the range
+[150,300]
+    ITS->SetRails(1);              // 1 --> rails in ; 0 --> rails out
+    ITS->SetCoolingFluid(1);       // 1 --> water ; 0 --> freon
 
-    AliITSvPPRsymm *ITS  = new AliITSvPPRsymm("ITS","New ITS PPR detailed version with symmetric services");
+    The first three functions are reserved to the ITS developpers and their
+values must not be modified at any time. The next four functions allow the user
+to  change the thickness of detectors and chips in the two SPD layers
+separately.  Detector thicknesses can go from 100 microns (TDR value) to 300
+microns (present default value). If a value outside this range is set a warning
+message is printed out and the default value (300 microns) is used. Chip
+thicknesses can go from 150 (TDR value) to 300 (present default value)
+microns. If a value outside this range is set a warning message is printed out
+and the default value (300 microns) is used. The last two function allow people
+to set the ITS rails in (=1) and out (=0) (the default is in) and the cooling
+fluid as water (=1, which is the default) or freon (=0). If the switch for
+rails is different from 0 and 1 a warning message is printed out and the
+default (rails in) is used. If the switch of the cooling fluid is different
+from 0 and 1 a warning message is printed out and the default (water) is used.
+
+     If you want to run the new PPR detailed geometry with symmetric services,
+the only uncommented lines must be:
+
+    AliITSvPPRsymm *ITS  = new AliITSvPPRsymm("ITS","New ITS PPR detailed version
+with symmetric services");
     ITS->SetMinorVersion(2);                                      
     ITS->SetReadDet(kFALSE);
     ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRsymm2.det");
-    ITS->SetThicknessDet1(300.);   // detector thickness on layer 1 must be in the range [100,300]
-    ITS->SetThicknessDet2(300.);   // detector thickness on layer 2 must be in the range [100,300]
-    ITS->SetThicknessChip1(300.);  // chip thickness on layer 1 must be in the range [150,300]
-    ITS->SetThicknessChip2(300.);  // chip thickness on layer 2 must be in the range [150,300]
+    ITS->SetThicknessDet1(300.);   // detector thickness on layer 1 must be in the
+range [100,300]
+    ITS->SetThicknessDet2(300.);   // detector thickness on layer 2 must be in the
+range [100,300]
+    ITS->SetThicknessChip1(300.);  // chip thickness on layer 1 must be in the range
+[150,300]
+    ITS->SetThicknessChip2(300.);  // chip thickness on layer 2 must be in the range
+[150,300]
     ITS->SetRails(1);              // 1 --> rails in ; 0 --> rails out
     ITS->SetCoolingFluid(1);       // 1 --> water ; 0 --> freon
 
+The functions are the same as in the asymmetric services' geometry. For their
+explanation and use see above.
+
 
 3. Simulation
 -------------  
 
-In order to run an event with a given ITS geometry put yourself in the directory
-macros and do the following:
+     In order to run an event with a given ITS geometry put yourself in the
+directory macros and do the following:
 
 - interactive run: start aliroot and the type the command "gAlice->Run()". 
   At the end of the run exit from aliroot with the command ".q".
@@ -164,78 +210,82 @@ macros and do the following:
   and error messages.
 
 In principle there is another way to run aliroot using the script alirun but
-this is not described here.
-Aliroot creates an output root file always called galice.root. If you want to
+this is not described here. Aliroot creates an output root file called
+galice.root (this can be changed from within the Config.C file). If you want to
 run the ITS reconstruction code, copy/move this file in the directory ITS and
-read the following instructions.  
-  
+read the following instructions.
 
 4. Cluster finding (fast)
 -------------------------
 
-Reconstructed points in the space can be calculated either smearing the Geant3 
+        Reconstructed points can be created quickly by smearing the Geant3 
 hits according with the various detector resolutions and applying the 
-thresholds for all detectors (fast reconstruction or, for short, "fast points"),
-or performing a cluster finding after the detector digitization (slow
-reconstruction or, for short, "slow points"). Fast point creation is described
-here while slow point one is described in sections 5 and 6.
-From now on, we assume that you are under the directory ITS. If it is the case
-you can do the following:
-
-- interactive run: start aliroot and the type the command 
-  ".x ITSHitsToFastPoints.C". At the end of the run exit from aliroot with the 
-  command ".q".
+thresholds for all detectors (fast reconstruction or, for short, "fast
+points"), or slowly, in more detail, by performing a cluster finding after the
+detector digitization (slow reconstruction or, for short, "slow points"). Fast
+point creation is described here while slow point one is described in sections
+5 and 6. From now on, we assume that you are under the directory ITS or that
+the $ALICE_ROOT/ITS is in your .rootrc MacroPath. If it is the case you can do
+the following, otherwise copy these macros to you local directory or prefix
+$ALICE_ROOT/ITS/ to the front of the macro names.
+
+- interactive run: start aliroot and type the command 
+  ".x ITSHitsToFastPoints.C".
+  At the end of the run exit from aliroot with the command ".q".
 - batch run: type the shell command 
-  "aliroot -q -b ITSHitsToFastPoints.C >& fileout &" where fileout is a name at 
-  your choice of the file where you want to store output and error messages. 
-  
-Fast points are written in the same galice.root file as you can see issuing the
-shell command "ls -l galice.root" before and after their creation and looking at
-the size of the root file. 
-By default, fast points are created for all kind of ITS subdetectors (SPD, SDD, 
-and SSD). This is done with the function call:
-
-  ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");  
+  "aliroot -q -b ITSHitsToFastPoints.C >& fileout &"
+  where fileout is a name at your choice of the file where you want to store
+  output and error messages.
+    Fast points are written in the same galice.root file as you can see issuing
+the shell command "ls -l galice.root" before and after their creation and
+looking at the size of the root file. 
+By default, fast points are created for all kind of ITS subdetectors (SPD, SDD,
+and SSD). This is done with the function call
+
+  ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
 
 in the macro ITSHitsToFastPoints.C. If you want to create fast points only for
 one type of subdectors you have to substitute the string "All" in the above 
 function call with "SPD", "SDD", or "SSD". Normal users are, however, strongly
-encouraged to create the "fast points" for all subdetectors at once not touching
-the macro ITSHitsToFastPoints.C.
-Fast points are intended only for tests. Normal users are also strongly 
-encouraged to run the complete ITS reconstruction described in the next two 
-sections. 
+encouraged to create the "fast points" for all subdetectors at once not
+touching the macro ITSHitsToFastPoints.C. Fast points are intended only for
+tests. Normal users are also strongly  encouraged to run the complete ITS
+reconstruction described in the next two sections. 
 
   
 5. Digitization
 ---------------
 
-In order to run the ITS digitization, put yourself under ITS and do the
-following:
+        The so called "Slow simulation" realy takes two parts. First digits are
+created (Digitization) and then these digits are read back in and are
+reconstructed to form RecPoints (Cluster finding). In order to run the ITS
+digitization, put yourself in ITS directory, or if you have the ITS directory
+in your .rootrc MacroPath, a working directory and do thefollowing:
 
-- interactive run: start aliroot and the type the command 
-  ".x ITSHitsToDigits.C". At the end of the run exit from aliroot with the 
-  command ".q".
+- interactive run: start aliroot and type the command 
+  ".x ITSHitsToDigits.C".
+  At the end of the run exit from aliroot with the command ".q".
 - batch run: type the shell command 
-  "aliroot -q -b ITSHitsToDigits.C >& fileout &" where fileout is a name at 
-  your choice of the file where you want to store output and error messages. 
-  
+  "aliroot -q -b ITSHitsToDigits.C >& fileout &" 
+  where fileout is a name at your choice of the file where you want to store
+  output and error messages.
+
 Digits are written in the same galice.root file as you can see issuing the
-shell command "ls -l galice.root" before and after their creation and looking at
-the size of the root file. 
-By default, digits are created for all kind of ITS subdetectors (SPD, SDD, 
-and SSD). This is done with the function call:
+shell command "ls -l galice.root" before and after their creation and looking
+at the size of the root file. By default, digits are created for all kind of
+ITS subdetectors (SPD, SDD, and SSD). This is done with the function call:
 
   ITS->HitsToDigits(nev,nbgr_ev,size," ","All"," "); 
 
 in the macro ITSHitsToDigits.C. If you want to create digits only for
 one type of subdectors you have to substitute the string "All" in the above 
 function call with "SPD", "SDD", or "SSD". Normal users are, however, strongly
-encouraged to run the ITS digitization for all subdetectors at once not touching
-the macro ITSHitsToDigits.C.
-By default the so-called "Dubna simulation" of the pixel detectors is 
-performed. In order to run the "Bari simulation" of the pixel detectors (waiting
-for the merging of the two) you have to use the macro ITSHitsToDigitsBari.C.
+encouraged to run the ITS digitization for all subdetectors at once not
+touching the macro ITSHitsToDigits.C. By default the so-called "Dubna
+simulation" of the pixel detectors is performed. In order to run the "Bari
+simulation" of the pixel detectors (waiting for the merging of the two) you
+have to use the macro ITSHitsToDigitsBari.C.
 
 
 6. Cluster finding (slow)
@@ -245,17 +295,17 @@ In order to perform the cluster finding and create the "slow points", do the
 following:
 
 - interactive run: start aliroot and the type the command 
-  ".x ITSDigitsToClusters.C". At the end of the run exit from aliroot with the 
-  command ".q".
+  ".x ITSDigitsToClusters.C".
+  At the end of the run exit from aliroot with the command ".q".
 - batch run: type the shell command 
-  "aliroot -q -b ITSDigitsToClusters.C >& fileout &" where fileout is a name at 
-  your choice of the file where you want to store output and error messages. 
-  
+  "aliroot -q -b ITSDigitsToClusters.C >& fileout &" 
+  where fileout is a name at  your choice of the file where you want to store
+  output and error messages.
+
 Slow points are written in the same galice.root file as you can see issuing the
-shell command "ls -l galice.root" before and after their creation and looking at
-the size of the root file. 
-By default, slow points are created for all kind of ITS subdetectors (SPD, SDD, 
-and SSD). This is done with the function call:
+shell command "ls -l galice.root" before and after their creation and looking
+at the size of the root file. By default, slow points are created for all kind
+of ITS subdetectors (SPD, SDD, and SSD). This is done with the function call:
 
   ITS->DigitsToRecPoints(nev,last_entry,"All"); 
 
@@ -263,71 +313,99 @@ in the macro ITSDigitsToClusters.C. If you want to create slow points only for
 one type of subdectors you have to substitute the string "All" in the above 
 function call with "SPD", "SDD", or "SSD". Normal users are, however, strongly
 encouraged to create the slow points for all subdetectors at once not touching
-the macro ITSDigitsToClusters.C.
-By default the so-called "Dubna reconstruction" of the pixel detectors is 
-performed. In order to run the "Bari reconstruction" of the pixel detectors 
-(waiting for the merging of the two) you have to use the macro 
-ITSDigitsToClustersBari.C.
+the macro ITSDigitsToClusters.C. By default the so-called "Dubna
+reconstruction" of the pixel detectors is performed. In order to run the "Bari
+reconstruction" of the pixel detectors (waiting for the merging of the two) you
+have to use the macro ITSDigitsToClustersBari.C.
 
 
 7. Useful test macros
 ---------------------
 
-Once you have created digits and slow points, you can now run some useful test
-macros. They can also be used as useful starting point to understand how to read
-hits, digits and points back from the galice.root file. The macros are listed 
-below together with some help on their use. Some of the macros produce output
-graphics, so the suggestion is to run them interactively starting aliroot and
-then typing the command ".x macro.C" at the aliroot prompt.
-
-- SPDclusterTest.C: this macro opens the galice.root file, reads the
-reconstructed points and plots the SPD resolution both for layer 1 and 2, in Z
-and Rphi direction. Moreover, it creates the Root file SPD_his.root which 
-contains some useful histograms and nt-uples which can be read back with the 
-macro SPD_ntuple.C.
-
-- ITSreadClustTestSPD.C: this macros opens the galice.root, reads SPD digits and
-prints them on the screen.
-
-- ITSsddanalysis.C: this macro opens the galice.root file, performs some
-analysis of the drift detectors and creates the Root output file 
-SDD_histos_test.root. This output file can then be read by the macro 
-ITSsddtest.C which create some PostScript file containing reference histograms.
-
-- SSDrecpointTest.C: this macro opens the galice.root file, reads the
-reconstructed points and plots the SSD resolution both for layer 5 and 6, in Z
-and Rphi direction. Moreover, it creates the Root file SSD_his.root which 
-contains some useful histograms and nt-uples which can be read back with the 
-macro SSD_ntuple.C.
-
-- ITSreadRecPointsTest.C: this macro opens the galice.root file, reads the
-reconstructed points and prints them on the screen for all ITS modules. For each
-reconstructed point the following quantities are printed: the recpoint index in
-the module (0, 1, 2, ...), the X coordinate in the local reference system of the
+        Once you have created digits and slow points, you can now run some
+useful test macros. They can also be used as useful starting point to
+understand how to read hits, digits and points back from the galice.root
+file. The macros are listed below together with some help on their use. Some of
+the macros produce graphical output, so the suggestion is to run them
+interactively from within aliroot by then typing a command like ".x macro.C" at
+the aliroot prompt.
+
+- new TBrower();
+  This isn't a macro, so you just type it at the aliroot prompt. From here you
+can open up root files and view some of what they contain. This is done by
+cliking on a root file, and then cliking on the different "folders" withing
+that root file. Data kept in TTrees can be histogramed symply by double cliking
+on the appropreate data member. Other parts of the root file can also be
+inspected in a simular way, right butting cliking for example.
+
+- SPDclusterTest.C:
+  This macro opens the galice.root file, reads the reconstructed points and
+plots the SPD resolution both for layer 1 and 2, in Z and Rphi
+direction. Moreover, it creates the Root file SPD_his.root which contains
+some useful histograms and nt-uples which can be read back with the macro
+SPD_ntuple.C. 
+
+- SPDclusterTestBari.C:
+  This macro opens the galice.root file, reads the reconstructed points and
+plots several distributions relative to the Bari simulation of the SPD. 
+Moreover, it creates the Root file SPD_his_bari.root which contains some useful
+histograms and nt-uples. This macro should also work with both SPD simulations.
+
+- ITSreadClustTestSPD.C:
+  This macros opens the galice.root, reads SPD digits and prints them on the 
+screen.
+
+- ITSsddanalysis.C:
+  This macro opens the galice.root file, performs some analysis of the 
+drift detectors and creates the Root output file SDD_histos_test.root. This 
+output file can then be read by the macro ITSsddtest.C which create some 
+PostScript file containing reference histograms.
+
+- SSDrecpointTest.C:
+  This macro opens the galice.root file, reads the reconstructed points and 
+plots the SSD resolution both for layer 5 and 6, in Z and Rphi direction. 
+Moreover, it creates the Root file SSD_his.root which contains some useful 
+histograms and nt-uples which can be read back with the macro SSD_ntuple.C.
+
+- ITSreadRecPointsTest.C:
+  This macro opens the galice.root file, reads the reconstructed points and 
+prints them on the screen for all ITS modules. For each reconstructed point 
+the following quantities are printed: the recpoint index in the module 
+(0, 1, 2, ...), the X coordinate in the local reference system of the
 module, the Z coordinate in the local reference system of the module, three
 integers indicating the id numbers of the tracks contributing to that
-reconstructed point. Only positive numbers, of course, are real tracks. Negative
-numbers (-1 for SPD, -3 for SDD, and -2 for SSD) are just there to fill the 
-vector of track id's associated to a given reconstructed point.
+reconstructed point. Only positive numbers, of course, are real tracks. 
+Negative numbers (-1 for SPD, -3 for SDD, and -2 for SSD) are just there to 
+fill the vector of track id's associated to a given reconstructed point.
 
-- ITSgeoplot.C: this macro opens the galice.root file, reads hits, digits and
-recpoint for SPD, SDD and SSD and plots them in the global reference system.
+- ITSgeoplot.C:
+  This macro opens the galice.root file, reads hits, digits and recpoint for 
+SPD, SDD and SSD and plots them in the global reference system.
 
-- ITSReadPlotData.C: this macros opens the galice.root file and then prompts the
-user for an ITS module identified by layer, ladder and detector or by the unique
-ID. That module is then searched for hits, digits and 
-reconstructed points and they are plotted in the local reference frame of the
-module with different symbols.
+- ITSReadPlotData.C:
+  This macros opens the galice.root file and then prompts the user for an ITS 
+module identified by layer, ladder and detector or by the unique ID. That 
+module is then searched for hits, digits and reconstructed points and they 
+are plotted in the local reference frame of the module with different symbols.
 
 
 8. Occupancy
 ------------
 
-ITS occupancy can be evaluated by the macro ITSOccupancy.C which opens the
-galice.root files, reads the number of digits for all subdetectors, and
-calculate the occupancy as the ratio of the "fired" digits over the total number
-of digits. Plots are also created showing the occupancy as a function of the Z
-coordinate along the beams' axis. In order to increase the speed of the 
+        ITS occupancy can be evaluated by the macro ITSOccupancy.C which opens
+the galice.root files, reads the number of digits for all subdetectors, and
+calculate the occupancy as the ratio of the "fired" digits over the total
+number of digits. Plots are also created showing the occupancy as a function of
+the Z coordinate along the beams' axis. In order to increase the speed of the 
 calculation the macro is compilable. See the comments in the code about how to 
 run it.
 
+
+9. Primary vertex
+-----------------
+
+        The ITS code contains a full 3D primary vertex finder. The relative
+class is AliITSVertex and an example of its use is contained in the test macro
+VertexMacro.C. The output of the primary vertex finder will eventually be 
+incorporated in the tracking algorithms.
+
diff --git a/ITS/SPDclusterTestBari.C b/ITS/SPDclusterTestBari.C
new file mode 100644 (file)
index 0000000..efaf903
--- /dev/null
@@ -0,0 +1,863 @@
+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();
+}
+
index 1306102..f3e9294 100644 (file)
@@ -142,9 +142,10 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
    //AliITSgeom *aliitsgeo = ITS->GetITSgeom();
    AliITSgeom *geom = ITS->GetITSgeom();
 
+
    //Int_t cp[8]={0,0,0,0,0,0,0,0};
 
-   cout << "SSD" << endl;
+   //cout << "SSD" << endl;
 
    AliITSDetType *iDetType=ITS->DetType(2);
    AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
@@ -155,13 +156,13 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
 
    TClonesArray *dig2  = ITS->DigitsAddress(2);
    TClonesArray *recp2  = ITS->ClustersAddress(2);
-   //   AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2,recp2);
    AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
    ITS->SetReconstructionModel(2,rec2);
+
    // test
-   printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz());
+   printf("SSD dimensions %f %f %f \n",seg2->Dx(),seg2->Dz(),seg2->Dy());
    printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
-
+   Float_t ylim = seg2->Dy()/2 - 12;
    
 //
 //   Loop over events
@@ -236,8 +237,8 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
      printf("det type %d first2, last2 %d %d \n",2,first2,last2);
 
      // module loop for the SSD
-     //for (mod=first2; mod<last2+1; mod++) {  // for the "ALL" option
-     for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option
+     for (mod=first2; mod<last2+1; mod++) {  // for the "ALL" option
+       //for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option
 
        TTree *TR = gAlice->TreeR();
        Int_t nentrec=TR->GetEntries();
@@ -255,8 +256,8 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
 
        Int_t nrecp = ITSrec->GetEntries();
        totpoints += nrecp;
-       if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
-       //if (!nrecp) continue;
+       //if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
+       if (!nrecp) continue;
        Int_t nclusters = ITSclu->GetEntries();
        totclust += nclusters;
        //if (nclusters) printf("Found %d clusters for module %d\n",nrecc,mod);
@@ -270,7 +271,7 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
        //       printf("Mod: %X\n",Mod);
        Int_t nhits = Mod->GetNhits();
        Float_t epart = 0;
-       //cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
+       cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
 
        // ---------------- cluster/hit analysis ---------------------
 
@@ -304,7 +305,6 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
         Int_t tr2 = itsPnt->GetLabel(ii);
          Int_t ii = 2;
         Int_t tr3 = itsPnt->GetLabel(ii);
-        //cout<<"pnt,ntrover,tr1,tr2,tr3 ="<<pnt<<","<<ntrover<<","<<tr1<<","<<tr2<<","<<tr3<<endl;
         // fill ntuple2
             ntuple2_st.nxP = nxP;
              ntuple2_st.nxN = nxN;
@@ -346,12 +346,11 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
    //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
    //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 
-          //cout<<"hit,hitstat,track,partcode,parent ="<<hit<<","<<hitstat<<","<<track<<","<<partcode<<","<<parent<<endl;
            Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
                                                      // vertex
           pmod *= 1.0e+3;
 
-         if(hitstat == 66 && yhit < -146.) {  // entering hit
+         if(hitstat == 66 && yhit < -ylim) {  // entering hit
            xhit0 = xhit;
            yhit0 = yhit;
            zhit0 = zhit;
@@ -369,7 +368,6 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
          // Consider the hits only with the track number equaled to one
          // of the recpoint
          if((track == tr1) || (track == tr2) || (track == tr3)) flagtrack = 1;
-         //cout<<"!!Ok: track,tr1,tr2,tr3 ="<<track<<","<<tr1<<","<<tr2<<","<<tr3<<endl;
 
          if(flagtrack == 1) {     // the hit corresponds to the recpoint
 
@@ -405,8 +403,6 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
            }
             pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
 
-            //cout<<"lay,pnt,hit,xmed,xrec,xdif,zmed,zrec,zdif ="<<hitlayer<<","<<pnt<<","<<hit<<","<<xmed<<","<<xrec<<","<<xdif<<","<<zmed<<","<<zrec<<","<<zdif<<endl;
-
             // fill ntuple
              ntuple_st.lay = hitlayer;
             ntuple_st.nxP = nxP;