ITS/ITShitOccupancy.C

   1 /*
   2 // Some time it is nice to compile the code to find all of the typos.
   3 #include <stdio.h>
   4 #include <math.h>
   5 #include "TROOT.h"
   6 #include "TSystem.h"
   7 #include "TClassTable.h"
   8 #include "TFile.h"
   9 #include "AliITS.h"
  10 #include "TTree.h"
  11 #include "TClonesArray.h"
  12 #include "TArray.h"
  13 #include "TCanvas.h"
  14 #include "AliRun.h"
  15 #include "AliITSgeom.h"
  16 #include "TParticlet.h"
  17
  18 extern TSystem     *gSystem;
  19 extern TROOT       *gROOT;
  20 extern TClassTable *gClassTable;
  21 extern AliRun      *gAlice;
  22 */
  23
  24 Double_t dEtadZ(Double_t z,Double_t r){
  25     Float_t a;
  26
  27     a = TMath::Sqrt(z*z+r*r);
  28     if(a!=0.0) a = 1./a;
  29     return(a);
  30 }
  31
  32 Double_t PdEtadZ(Double_t *z,Double_t *r){
  33     return(dEtadZ(*z,*r));
  34 }
  35
  36 Double_t ZfromEtaR(Double_t eta,Double_t r){
  37     Double_t z;
  38
  39     z   = 2.0*TMath::ATan(TMath::Exp(-eta));
  40     z   = r/TMath::Tan(z);
  41     return(z);
  42 }
  43
  44 Double_t EtafromZr(Double_t z,Double_t r){
  45     Double_t a,eta;
  46
  47     a   = 0.5*TMath::ATan2(r,z);    // Tan2(y,x) = y/x
  48     eta = -TMath::Log(TMath::Abs(TMath::Tan(a)));
  49     if(a<0.0) eta = -eta;
  50     return (eta);
  51 }
  52
  53 Double_t PEtafromZr(Double_t *z,Double_t *r){
  54     return (EtafromZr(*z,*r));
  55 }
  56
  57 Double_t Weight(Double_t *a,Double_t *b){
  58     Double_t eta,z,Etap,rt,dEta,dEtap,weight,REtap;
  59
  60     eta    = a[0];
  61     Etap   = a[1];
  62     rt     = b[0];
  63     dEta   = b[1];
  64     dEtap  = b[2];
  65     REtap  = b[3];
  66     z      = ZfromEtaR(eta,rt);
  67     weight = TMath::Abs(dEtadZ(z,rt));
  68     weight = weight*REtap/(rt*2.0*TMath::Pi()*dEta*dEtap);
  69     return weight;
  70 }
  71
  72 void ITShitOccupancy (const char *filename="galice.root",Int_t evNumber=0){
  73 /////////////////////////////////////////////////////////////////////////
  74 //
  75 //   Written by Roberto Barbera
  76 //   Modified by Bjorn S. Nilsen May 19 1999
  77 //
  78 //   This macro is a small example of a ROOT macro
  79 //   illustrating how to read the output of aliroot
  80 //   and fill some histograms. The Occupancy values produce are only
  81 //   approximate and are known to have problems. No warranty of any kind
  82 //   may be give or assumed.
  83 //
  84 //     Root > .L ITShitOccupancy.C//this loads the macro in memory
  85 //     Root > ITShitOccupancy();  //by default process first event
  86 //     Root > ITShitOccupancy("galiceA.root");  //process file galiceA.root
  87 //     Root > ITShitOccupancy("galiceA.root",2);//process file galiceA.root
  88 //                                                 third event
  89 //
  90 //     Root > .x ITShitOccupancy.C(); //by default process first event
  91 //     Root > .x ITShitOccupancy.C("galiceA.root"); //process file galiceA.root
  92 //     Root > .x ITShitOccupancy.C("galiceA.root",2);//process file
  93 //                                                     galiceA.root third event
  94 //     The test figures for this macro are kept in
  95 // $ALICE_ROOT/ITS/figures/ITShitOccupancy.figures as a series of post-
  96 // script files and one ITShitOccupancy.root file. This is because ther
  97 // are so many figures produced.
  98 //
  99 ////////////////////////////////////////////////////////////////////////
 100 //
 101 // The code computes the Occupancy of HITS in the different ITS detectors
 102 // as a function of z. This does not represent the true occupancy expected
 103 // in the ITS since there will be more than one hit per charge cluster in
 104 // the different detector of the ITS. To do a proper occupancy calculation
 105 // a full simulation of all of the detectors in the ITS will be needed. In
 106 // addition a proper background needs to be used.
 107 //
 108 //  This code has been modified to be as compatible with the standard (May 26
 109 //  1999) version of the official released ITS code as possible. Functions
 110 //  introduce in the unofficial release of the ITS code (March 199) have
 111 //  been commented out and replaced by the variables defined in the ITShits
 112 //  class. As soon as these functions are officially part of the aliroot
 113 //  release. These functions should be used instead to maximise future
 114 //  portability.
 115 //
 116 //  As with all ITS software produced under the deadline of the ITS TDR, there
 117 //  is no real guarantee that this code will function, produce meaning full
 118 //  results, solve in any way any possible problem or question, or be supported
 119 //  in any kind. Questions or comments can be addressed to
 120 //  Nilsen@mps.ohio-state.edu (for as long as I am there) or
 121 //  Roberto.Barbera@ct.infn.it. A response may never come.
 122 //
 123 //  This exclamation should be copy written and thereby not used my any other
 124 //  software writer, especially those writing code for ALICE. We expect
 125 //  everyone else's code to work perfectly and do the job as advertised, either
 126 //  explicitly, by word of mouth, or by any other mystical means.
 127 //
 128 ////////////////////////////////////////////////////////////////////////
 129 // Dynamically link some shared libs
 130    if (gClassTable->GetID("AliRun") < 0) {
 131       gROOT->LoadMacro("loadlibs.C");
 132       loadlibs();
 133    } // end if gClassTable...
 134 //
 135 // Connect the Root Galice file containing Geometry, Kine and Hits
 136    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
 137    if (!file) file = new TFile(filename);
 138    printf("Reading from root file %s\n",filename);
 139 //
 140 // Get AliRun object from file or create it if not on file
 141    if (!gAlice) {
 142       gAlice = (AliRun*)file->Get("gAlice");
 143       if (gAlice) printf("AliRun object found on file\n");
 144       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 145    } // end if !gAlice
 146 //
 147 // Get pointers to Alice detectors and Hits containers
 148    Int_t  nparticles  = gAlice->GetEvent(evNumber);
 149    if (nparticles <= 0) return;
 150    AliITS       *ITS       = gAlice->GetDetector("ITS");
 151    if(!ITS) {printf("ITS detector not found. Exiting\n"); return;}
 152    TClonesArray *Particles = gAlice->Particles();
 153    TClonesArray *ITShits   = ITS->Hits();
 154    TTree        *TH        = gAlice->TreeH();
 155    Int_t        ntracks    = TH->GetEntries();
 156    AliITSgeom   *gm        = ITS->GetITSgeom();
 157    Int_t        Nlayers    = gm->GetNlayers();
 158 //
 159    printf("%d primary tracks and their secondary tracks read in\n",ntracks);
 160 //
 161 // Import the Kine and Hits Trees for the event evNumber in the file
 162    Float_t x,y,z,mass,e,r,rt;
 163    Float_t px,py,pz,pt;
 164    Float_t destep;
 165    Float_t theta,phi;
 166    Float_t eta,etap;
 167    Float_t dzeta,weight;
 168    Int_t idpart;
 169 //
 170    Float_t ztot[]  = {33.52,// total length of layer 1
 171                       33.52,// total length of layer 2
 172                       46.10,// total length of layer 3
 173                       60.77,// total length of layer 4
 174                       90.42,// total length of layer 5
 175                      101.95};// total length of layer 6
 176    Float_t raprx[] = {4.0,7.0,14.9,23.8,39.1,43.6}; // aproximate layer radii
 177 //
 178    Float_t totdet[]={256*256*4*20,//total number of channels layer 1 =  5242880
 179                      256*256*4*40,//total number of channels layer 2 = 10485760
 180                    2*256*256*6*14,//total number of channels layer 3 = 11010048
 181                    2*256*256*8*22,//total number of channels layer 4 = 23068672
 182                         768*23*34,//total number of channels layer 5 =  1201152
 183                         768*26*38};//total number of channels layer 6 = 1517568
 184 //
 185    Float_t clusize,clusize_z,clusize_rphi;
 186 //
 187    Float_t threshold[]={8.0e-6,8.0e-6, // noise= 200 e- ; thr = 2000 e- ;
 188                         3.0e-6,3.0e-6, // noise= 250 e- ; thr= 3 sigma noise
 189                         6.0e-6,6.0e-6};// noise= 500 e- ; thr= 3 sigma noise
 190 //
 191    Int_t     nhitlayer[6];
 192    Int_t     ntrklayer[6];
 193    Int_t     layer,ladder,det;
 194    Int_t     nbytes = 0;
 195    Int_t     i,j,jj,hit,ipart,ipartold,iparent;
 196    Int_t     nhits;
 197    Int_t     sector,plane;
 198    TParticle *particle,*pparticle;
 199    AliITShit *itsHit;
 200
 201 // Create histograms
 202
 203    Float_t etamin = -1.;
 204    Float_t etamax = 1.;
 205    Float_t ptmin  = 0.;
 206    Float_t ptmax  = 3.;
 207    Int_t   nbin   = 40;
 208    Int_t   nbinz  = 200;
 209    Float_t deta   = (etamax-etamin)/nbin;
 210    Float_t EtapMax=2.436246,EtapMin=-2.436246; // from Config.C file
 211    TVector etabin(nbin+1);
 212    TVector zbin(nbin+1);
 213
 214 // Open a root file to save histograms in
 215    char *Hfilename = "ITShitOccupancy.root";
 216    TFile *Hfile    = (TFile *)gROOT->GetListOfFiles()->FindObject(Hfilename);
 217    if(Hfile) Hfile->Close();
 218    Hfile = new TFile(Hfilename,"RECREATE","Histograms for ITS Occupancy");
 219    printf("Writting histograms to root file %s\n",Hfilename);
 220 //
 221    TH2F *hzetap[6];
 222    TH1D *hzetap_z[6],*hzetap_etap[6];
 223    TH2F *hzetapden[6];
 224    TH1D *hzetapden_z[6],*hzetapden_etap[6];
 225    TH2F *hEtaEtapden[6];
 226    TH1D *hEtaEtapden_Eta[6],*hEtaEtapden_Etap[6];
 227    TH2F *hEtaEtapdenWeight[6];
 228    TH1D *hEtaEtapdenWeight_Eta[6];
 229    TF2  *fweight;
 230    TF1  *fdEtadZ,*fEta;
 231    TH1F *hdEtadZ[6],*hEta[6];
 232
 233    { // "Book" histograms
 234        char histid[7],htit[40];
 235        for(i=0;i<Nlayers;i++){
 236            sprintf(histid,"hzetap%1.1d",i+1);
 237            sprintf(htit,"Hits for Layer %1.1d",i+1);
 238            hzetap[i] = new TH2F(histid,htit,nbinz,-0.6*ztot[i],0.6*ztot[i],
 239                                 20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
 240            hzetap[i]->SetXTitle("Z (cm)");
 241            hzetap[i]->SetYTitle("Beam particle pseudorapidity");
 242            hzetap[i]->Sumw2();
 243            //
 244            sprintf(histid,"hzetapden%1.1d",i+1);
 245            sprintf(htit,"Density of Hits for Layer %1.1d",i+1);
 246            hzetapden[i] = new TH2F(histid,htit,nbinz,-0.6*ztot[i],0.6*ztot[i],
 247                                20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
 248            hzetapden[i]->Sumw2();
 249            hzetapden[i]->SetXTitle("Z (cm)");
 250            hzetapden[i]->SetYTitle("Beam particle pseudorapidity");
 251            nhitlayer[i] = 0.0;
 252            ntrklayer[i] = 0.0;
 253            //
 254            sprintf(histid,"hEtaEtapden%1.1d",i+1);
 255            sprintf(htit,"Hits for Layer %1.1d by Eta",i+1);
 256            hEtaEtapden[i] = new TH2F(histid,htit,nbinz,
 257                                   -EtafromZr(0.6*ztot[i],raprx[i]),
 258                                    EtafromZr(0.6*ztot[i],raprx[i]),
 259                                 20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
 260            hEtaEtapden[i]->SetXTitle("Track Pseudorapidity");
 261            hEtaEtapden[i]->SetYTitle("Beam particle pseudorapidity");
 262            hEtaEtapden[i]->Sumw2();
 263        } // end for i
 264
 265        sprintf(histid,"fweight");
 266        fweight = new TF2(histid,Weight,-EtafromZr(0.6*ztot[0],raprx[0]),
 267                                    EtafromZr(0.6*ztot[0],raprx[0]),
 268                          floor(EtapMin)-1.0,ceil(EtapMax)+1.0,4);
 269        Double_t params[4];
 270        sprintf(histid,"fdEtadZ");
 271        fdEtadZ = new TF1(histid,PdEtadZ,floor(EtapMin)-1.0,
 272                          ceil(EtapMax)+1.0,1);
 273        sprintf(histid,"fEta");
 274        fEta = new TF1(histid,PEtafromZr,floor(EtapMin)-1.0,
 275                       ceil(EtapMax)+1.0,1);
 276        for(i=0;i<Nlayers;i++){
 277            params[0] = raprx[i];
 278            params[1] = hEtaEtapden[i]->GetXaxis()->GetBinWidth(1);
 279            params[2] = hEtaEtapden[i]->GetYaxis()->GetBinWidth(1);
 280            params[3] = EtapMax-EtapMin;
 281            fweight->SetParameters(params);
 282            sprintf(histid,"hEtaEtapdenWeight%1.1d",i+1);
 283            sprintf(htit,"Weight for Layer %1.1d by Eta",i+1);
 284            hEtaEtapdenWeight[i] = new TH2F(histid,htit,nbinz,
 285                                   -EtafromZr(0.6*ztot[i],raprx[i]),
 286                                    EtafromZr(0.6*ztot[i],raprx[i]),
 287                                 20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
 288            hEtaEtapdenWeight[i]->SetXTitle("Track Pseudorapidity");
 289            hEtaEtapdenWeight[i]->SetYTitle("Beam particle pseudorapidity");
 290            hEtaEtapdenWeight[i]->Sumw2();
 291            hEtaEtapdenWeight[i]->Eval(fweight);
 292 //
 293            fdEtadZ->SetParameters(params);
 294            sprintf(histid,"hdEtadZ%1.1d",i+1);
 295            sprintf(htit,"dEtadZ(z) for Layer %1.1d",i+1);
 296            hdEtadZ[i] = new TH1F(histid,htit,nbinz,
 297                                   -0.6*ztot[i],0.6*ztot[i]);
 298            hdEtadZ[i]->SetXTitle("Z (cm)");
 299            hdEtadZ[i]->Sumw2();
 300            hdEtadZ[i]->Eval(fdEtadZ);
 301 //
 302            fEta->SetParameters(params);
 303            sprintf(histid,"hEta%1.1d",i+1);
 304            sprintf(htit,"Eta(z) for Layer %1.1d",i+1);
 305            hEta[i] = new TH1F(histid,htit,nbinz,
 306                                   -0.6*ztot[i],0.6*ztot[i]);
 307            hEta[i]->SetXTitle("Z (cm)");
 308            hEta[i]->Sumw2();
 309            hEta[i]->Eval(fEta);
 310        } // end for i
 311
 312    } //  end "Book" histograms
 313 // Start loop on tracks in the hits container
 314
 315    Double_t xa[2],par[4];
 316
 317    for (Int_t track=0; track<ntracks;track++) {
 318      gAlice->ResetHits();
 319      nbytes += TH->GetEvent(track);
 320      nhits   = ITShits->GetEntriesFast();
 321      for (hit=0;hit<nhits;hit++) {
 322          itsHit   = (AliITShit*)ITShits->UncheckedAt(hit);
 323          destep   = itsHit->GetIonization();
 324          // With this new version, to be able to do proper detector
 325          // simulations, the requirment that a "hit" leave some
 326          // energy deposited has been removed.
 327          if(destep<=0.0) continue; // skip hits without energy loss.
 328          ipart    = itsHit->fTrack;
 329          particle = (TParticle*)Particles->UncheckedAt(ipart);
 330          iparent  = particle->GetFirstMother();
 331          pparticle= (TParticle*)Particles->UncheckedAt(ntracks-track);
 332          etap     = TMath::ATan2(pparticle->Pt(),pparticle->Pz());
 333          etap     = -TMath::Log(TMath::Tan(0.5*etap));
 334          itsHit->GetPositionG(x,y,z);
 335          x -= pparticle->Vx();  // correct for primary vertex position
 336          y -= pparticle->Vy();  // correct for primary vertex position
 337          z -= pparticle->Vz();  // correct for primary vertex position
 338          itsHit->GetMomentumG(px,py,pz);
 339          itsHit->GetDetectorID(layer,ladder,det);
 340          r        = sqrt(x*x+y*y+z*z);
 341          rt       = sqrt(x*x+y*y);
 342          theta    = TMath::ACos(z/r)*180./TMath::Pi();
 343          phi      = TMath::ATan2(y,x); if(phi<0.0) phi += TMath::Pi();
 344          eta      = EtafromZr(z,rt);
 345          pt       = TMath::Sqrt(px*px+py*py);
 346          if(ipart!=ipartold)ntrklayer[layer-1] += 1;
 347          ipartold = ipart;
 348          nhitlayer[layer-1] += 1;
 349
 350          switch (layer) {
 351         // Calculate average cluster sizes
 352          case 1:
 353              clusize_rphi = 1.4;
 354              if(TMath::Abs(eta)>=0.0 && TMath::Abs(eta)< 0.2) clusize_z=1.2;
 355              if(TMath::Abs(eta)>=0.2 && TMath::Abs(eta)< 0.4) clusize_z=1.2;
 356              if(TMath::Abs(eta)>=0.4 && TMath::Abs(eta)< 0.6) clusize_z=1.3;
 357              if(TMath::Abs(eta)>=0.6 && TMath::Abs(eta)< 0.8) clusize_z=1.4;
 358              if(TMath::Abs(eta)>=0.8 && TMath::Abs(eta)< 1.0) clusize_z=1.5;
 359              clusize = clusize_z*clusize_rphi;
 360          case 2:
 361              clusize_rphi = 1.4;
 362              if(TMath::Abs(eta)>=0.0 && TMath::Abs(eta)< 0.2) clusize_z=1.2;
 363              if(TMath::Abs(eta)>=0.2 && TMath::Abs(eta)< 0.4) clusize_z=1.2;
 364              if(TMath::Abs(eta)>=0.4 && TMath::Abs(eta)< 0.6) clusize_z=1.3;
 365              if(TMath::Abs(eta)>=0.6 && TMath::Abs(eta)< 0.8) clusize_z=1.4;
 366              if(TMath::Abs(eta)>=0.8 && TMath::Abs(eta)< 1.0) clusize_z=1.5;
 367              clusize = clusize_z*clusize_rphi;
 368          case 3:
 369              clusize_rphi = 8.0;
 370              clusize_z    = 3.0;
 371              clusize      = clusize_z*clusize_rphi;
 372          case 4:
 373              clusize_rphi = 8.0;
 374              clusize_z    = 3.0;
 375              clusize      = clusize_z*clusize_rphi;
 376          case 5:
 377              clusize_z = 1.0;
 378              if(pt>=0.0 && pt< 0.2) clusize_rphi=2.0;
 379              if(pt>=0.2 && pt< 0.4) clusize_rphi=1.9;
 380              if(pt>=0.4 && pt< 0.6) clusize_rphi=1.6;
 381              if(pt>=0.6 && pt< 0.8) clusize_rphi=1.6;
 382              if(pt>=0.8 && pt< 1.0) clusize_rphi=1.7;
 383              if(pt>=1.0 && pt< 1.5) clusize_rphi=1.7;
 384              if(pt>=1.5 && pt< 2.5) clusize_rphi=1.7;
 385              if(pt>=2.5           ) clusize_rphi=1.7;
 386              clusize = clusize_z*clusize_rphi;
 387          case 6:
 388              clusize_z = 1.0;
 389              if(pt>=0.0 && pt< 0.2) clusize_rphi=2.0;
 390              if(pt>=0.2 && pt< 0.4) clusize_rphi=1.9;
 391              if(pt>=0.4 && pt< 0.6) clusize_rphi=1.6;
 392              if(pt>=0.6 && pt< 0.8) clusize_rphi=1.6;
 393              if(pt>=0.8 && pt< 1.0) clusize_rphi=1.7;
 394              if(pt>=1.0 && pt< 1.5) clusize_rphi=1.7;
 395              if(pt>=1.5 && pt< 2.5) clusize_rphi=1.7;
 396              if(pt>=2.5           ) clusize_rphi=1.7;
 397              clusize = clusize_z*clusize_rphi;
 398          } // end switch layer
 399 //       weightE = clusize*ztot[layer-1]/totdet[layer-1];
 400          xa[0]  = eta; xa[1]  = etap;
 401          par[0] = rt;
 402          par[1] = hEtaEtapden[layer-1]->GetXaxis()->GetBinWidth(1);
 403          par[2] = hEtaEtapden[layer-1]->GetYaxis()->GetBinWidth(1);
 404          par[3] = EtapMax-EtapMin;
 405          weight = Weight(xa,par);
 406          hEtaEtapden[layer-1]->Fill(eta,etap,weight);
 407          weight = 1.;
 408          hzetap[layer-1]->Fill(z,etap,weight);
 409          weight = par[3]/(rt*2.0*TMath::Pi()*
 410                        hzetapden[layer-1]->GetXaxis()->GetBinWidth(1)*
 411                        hzetapden[layer-1]->GetYaxis()->GetBinWidth(1));
 412          hzetapden[layer-1]->Fill(z,etap,weight);
 413        } // end for hit
 414    } // end for track
 415
 416    for(i=0;i<Nlayers;i++)printf("No. of hits   in layer %d = %d\n",i+1,nhitlayer[i]);
 417    for(i=0;i<Nlayers;i++)printf("No. of tracks in layer %d = %d\n",i+1,ntrklayer[i]);
 418
 419 //Create canvases, set the view ranges and show histograms
 420
 421      TCanvas *canvas = new TCanvas("canvas","ITShitOccupancy",1);
 422
 423      Bool_t printit = kTRUE;
 424      char psfilename[40];
 425      for(i=0;i<Nlayers;i++){
 426          hzetap[i]->Draw("colz");
 427          sprintf(psfilename,"ITShitOccupancy_%1.1d_z_etap.ps",i+1);
 428          if(printit) canvas->Print(psfilename);
 429          sprintf(psfilename,"hzetap_z%1.1d",i+1);
 430          hzetap_z[i]    = hzetap[i]->ProjectionX(psfilename,0,
 431                                                  hzetap[i]->GetNbinsY()+1,"E");
 432          hzetap_z[i]->SetXTitle("Z (cm)");
 433          hzetap_z[i]->SetYTitle("Number of Hits per cm");
 434          sprintf(psfilename,"hzetap_etap%1.1d",i+1);
 435          hzetap_etap[i] = hzetap[i]->ProjectionY(psfilename,0,
 436                                                  hzetap[i]->GetNbinsX()+1,"E");
 437          hzetap_etap[i]->SetXTitle("Beam particle pseudorapidity");
 438          hzetap_etap[i]->SetYTitle("Number of Hits");
 439          hzetap_z[i]->Draw();
 440          sprintf(psfilename,"ITShitOccupancy_%1.1d_z.ps",i+1);
 441          if(printit) canvas->Print(psfilename);
 442          hzetap_etap[i]->Draw();
 443          sprintf(psfilename,"ITShitOccupancy_%1.1d_etap.ps",i+1);
 444          if(printit) canvas->Print(psfilename);
 445 //
 446          hzetapden[i]->Draw("colz");
 447          sprintf(psfilename,"ITSHitDen_%1.1d_z_etap.ps",i+1);
 448          if(printit) canvas->Print(psfilename);
 449          sprintf(psfilename,"hzetapden_z%1.1d",i+1);
 450          hzetapden_z[i]    = hzetapden[i]->ProjectionX(psfilename,0,
 451                                         hzetapden[i]->GetNbinsY()+1,"E");
 452          hzetapden_z[i]->SetXTitle("Z (cm)");
 453          hzetapden_z[i]->SetYTitle("Number of Hits (cm^-2)");
 454          sprintf(psfilename,"hzetapden_etap%1.1d",i+1);
 455          hzetapden_etap[i] = hzetapden[i]->ProjectionY(psfilename,0,
 456                                       hzetapden[i]->GetNbinsX()+1,"E");
 457          hzetapden_etap[i]->SetXTitle("Beam particle pseudorapidity");
 458          hzetapden_etap[i]->SetYTitle("Number of Hits (cm^-2)");
 459          hzetapden_z[i]->Draw();
 460          sprintf(psfilename,"ITSHitDen_%1.1d_z.ps",i+1);
 461          if(printit) canvas->Print(psfilename);
 462          hzetapden_etap[i]->Draw();
 463          sprintf(psfilename,"ITSHitDen_%1.1d_etap.ps",i+1);
 464          if(printit) canvas->Print(psfilename);
 465 //
 466          hEtaEtapden[i]->Draw("colz");
 467          sprintf(psfilename,"ITShitDensity_%1.1d_Eta_Etap.ps",i+1);
 468          if(printit) canvas->Print(psfilename);
 469          sprintf(psfilename,"hEtaEtapden_Eta%1.1d",i+1);
 470          hEtaEtapden_Eta[i]    = hEtaEtapden[i]->ProjectionX(psfilename,0,
 471                                          hEtaEtapden[i]->GetNbinsY()+1,"E");
 472          hEtaEtapden_Eta[i]->SetXTitle("Pseudorapidity");
 473          hEtaEtapden_Eta[i]->SetYTitle("Number of Hits (cm^-2)");
 474          sprintf(psfilename,"hEtaEtapden_Etap%1.1d",i+1);
 475          hEtaEtapden_Etap[i] = hEtaEtapden[i]->ProjectionY(psfilename,0,
 476                                          hEtaEtapden[i]->GetNbinsX()+1,"E");
 477          hEtaEtapden_Etap[i]->SetXTitle("Beam particle pseudorapidity");
 478          hEtaEtapden_Etap[i]->SetYTitle("Number of Hits (cm^-2)");
 479          hEtaEtapden_Eta[i]->Draw();
 480          sprintf(psfilename,"ITShitDensity_%1.1d_Eta.ps",i+1);
 481          if(printit) canvas->Print(psfilename);
 482          hEtaEtapden_Etap[i]->Draw();
 483          sprintf(psfilename,"ITShitDensity_%1.1d_Etap.ps",i+1);
 484          if(printit) canvas->Print(psfilename);
 485 //
 486          hEtaEtapdenWeight[i]->Draw("colz");
 487          sprintf(psfilename,"ITShitDensityWeight_%1.1d_Eta_Etap.ps",i+1);
 488          if(printit) canvas->Print(psfilename);
 489          sprintf(psfilename,"hEtaEtapdenWeight_Eta%1.1d",i+1);
 490          hEtaEtapdenWeight_Eta[i] =
 491                                hEtaEtapdenWeight[i]->ProjectionX(psfilename,
 492                              0,hEtaEtapdenWeight[i]->GetNbinsY()+1,"E");
 493          hEtaEtapdenWeight_Eta[i]->SetXTitle("Pseudorapidity");
 494          hEtaEtapdenWeight_Eta[i]->SetYTitle("Weight (cm^-2)");
 495          hEtaEtapdenWeight_Eta[i]->Draw();
 496          sprintf(psfilename,"ITShitDensityWeight_%1.1d_Eta.ps",i+1);
 497          if(printit) canvas->Print(psfilename);
 498      } // end for i
 499 //
 500      Hfile->Write();
 501 //
 502      return;
 503 }