ITS/SPDclusterTestDubna.C

   1 #include "iostream.h"
   2
   3 void SPDclusterTestDubna (Int_t evNumber1=0,Int_t evNumber2=0)
   4 {
   5 /////////////////////////////////////////////////////////////////////////
   6 //   This macro is a small example of a ROOT macro
   7 //   illustrating how to read the output of GALICE
   8 //   and do some analysis.
   9 //
  10 /////////////////////////////////////////////////////////////////////////
  11
  12 // Dynamically link some shared libs
  13
  14    if (gClassTable->GetID("AliRun") < 0) {
  15       gROOT->LoadMacro("loadlibs.C");
  16       loadlibs();
  17    } else {
  18       delete gAlice;
  19       gAlice=0;
  20    }
  21
  22 // Connect the Root Galice file containing Geometry, Kine and Hits
  23
  24    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  25    if (!file) file = new TFile("galice.root");
  26    file->ls();
  27
  28 // Get AliRun object from file or create it if not on file
  29
  30    if (!gAlice) {
  31       gAlice = (AliRun*)file->Get("gAlice");
  32       if (gAlice) printf("AliRun object found on file\n");
  33       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  34    }
  35
  36           // ------------ Cluster and point analysis histogramms ------------
  37
  38           TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.);
  39           TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.);
  40           TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.);
  41           TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.);
  42           TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.);
  43           TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.);
  44           TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.);
  45           TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.);
  46
  47           TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.);
  48           TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.);
  49           TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.);
  50           TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.);
  51
  52           TH1F *Ptot1 = new TH1F("Ptot1","Total momentum (GeV/C) for layers 1",100,0.,5.);
  53           TH1F *Pz1 = new TH1F("Pz1","Pz (GeV/C) for layers 1",100,-5.,5.);
  54           TH1F *Theta1 = new TH1F("Theta1","Theta angle (rad) for layers 1",100,0.,4.);
  55           TH1F *Y1 = new TH1F("Y1","Rapidity for layers 1",100,-4.,4.);
  56           TH1F *Eta1 = new TH1F("Eta1","PseudoRapidity for layers 1",100,-4.,4.);
  57           TH1F *Y1Den = new TH1F("Y1Den","Rapidity for layers 1",100,-0.5,0.5);
  58           TH1F *Eta1Den = new TH1F("Eta1Den","PseudoRapidity for layers 1",100,-0.5,0.5);
  59           TH1F *Y1DenA = new TH1F("Y1DenA","Rapidity for layers 1",100,-0.5,0.5);
  60           TH1F *Eta1DenA = new TH1F("Eta1DenA","PseudoRapidity for layers 1",100,-0.5,0.5);
  61           TH1F *Phi1 = new TH1F("Phi1","Phi angle (rad) for layers 1",100,0.,7.);
  62           TH1F *Ptot2 = new TH1F("Ptot2","Total momentum (GeV/C) for layers 2",100,0.,5.);
  63           TH1F *Pz2 = new TH1F("Pz2","Pz (GeV/C) for layers 2",100,-5.,5.);
  64           TH1F *Theta2 = new TH1F("Theta2","Theta angle (rad) for layers 2",100,0.,4.);
  65           TH1F *Y2 = new TH1F("Y2","Rapidity for layers 2",100,-4.,4.);
  66           TH1F *Eta2 = new TH1F("Eta2","PseudoRapidity for layers 2",100,-4.,4.);
  67           TH1F *Y2Den = new TH1F("Y2Den","Rapidity for layers 2",100,-0.5,0.5);
  68           TH1F *Eta2Den = new TH1F("Eta2Den","PseudoRapidity for layers 2",100,-0.5,0.5);
  69           TH1F *Y2DenA = new TH1F("Y2DenA","Rapidity for layers 2",100,-0.5,0.5);
  70           TH1F *Eta2DenA = new TH1F("Eta2DenA","PseudoRapidity for layers 2",100,-0.5,0.5);
  71           TH1F *Phi2 = new TH1F("Phi2","Phi angle (rad) for layers 2",100,0.,7.);
  72
  73           // -------------- Create ntuples --------------------
  74           //  ntuple structures:
  75
  76           struct {
  77             Int_t lay;
  78             Int_t nx;
  79             Int_t nz;
  80             Int_t hitprim;
  81             Int_t partcode;
  82             Int_t ntrover;
  83             Float_t dx;
  84             Float_t dz;
  85             Float_t pmod;
  86           } ntuple_st;
  87
  88           struct {
  89             Int_t lay;
  90             Int_t lad;
  91             Int_t det;
  92             Int_t nx;
  93             Int_t nz;
  94             Int_t ntrover;
  95             Int_t noverlaps;
  96             Int_t noverprim;
  97             Float_t qcl;
  98             Float_t dx;
  99             Float_t dz;
 100             Float_t x;
 101             Float_t z;
 102           } ntuple1_st;
 103
 104           struct {
 105             //      Int_t lay;
 106             Int_t lay;
 107             Int_t nx;
 108             Int_t nz;
 109             Float_t x;
 110             Float_t z;
 111             Float_t qcl;
 112           } ntuple2_st;
 113
 114           ntuple = new TTree("ntuple","Demo ntuple");
 115           ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
 116           ntuple->Branch("nx",&ntuple_st.nx,"nx/I");
 117           ntuple->Branch("nz",&ntuple_st.nz,"nz/I");
 118           ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
 119           ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
 120           ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
 121           ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
 122           ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
 123           ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
 124
 125           ntuple1 = new TTree("ntuple1","Demo ntuple1");
 126           ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
 127           ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
 128           ntuple1->Branch("det",&ntuple1_st.det,"det/I");
 129           ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I");
 130           ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I");
 131           ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F");
 132           ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
 133           ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
 134           ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
 135           ntuple1->Branch("x",&ntuple1_st.x,"x/F");
 136           ntuple1->Branch("z",&ntuple1_st.z,"z/F");
 137           ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
 138           ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
 139
 140
 141           ntuple2 = new TTree("ntuple2","Demo ntuple2");
 142           //      ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
 143           ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
 144           ntuple2->Branch("x",&ntuple2_st.x,"x/F");
 145           ntuple2->Branch("z",&ntuple2_st.z,"z/F");
 146           ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I");
 147           ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I");
 148           ntuple2->Branch("qcl",&ntuple2_st.qcl,"qcl/F");
 149
 150 // ------------------------------------------------------------------------
 151 //
 152 //   Loop over events
 153 //
 154    for (int nev=0; nev<= evNumber2; nev++) {
 155      Int_t nparticles = gAlice->GetEvent(nev);
 156      cout << "nev         " << nev <<endl;
 157      cout << "nparticles  " << nparticles <<endl;
 158      if (nev < evNumber1) continue;
 159      if (nparticles <= 0) return;
 160
 161      TTree *TH = gAlice->TreeH();
 162      Int_t ntracks = TH->GetEntries();
 163      cout<<"ntracks "<<ntracks<<endl;
 164
 165 // Get pointers to Alice detectors and Digits containers
 166    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 167    TClonesArray *Particles = gAlice->Particles();
 168
 169    if (ITS) {
 170      // fill modules with sorted by module hits
 171      Int_t nmodules;
 172      ITS->InitModules(-1,nmodules);
 173      //    ITS->FillModules(nev,-1,evNumber2,nmodules," "," ");
 174      ITS->FillModules(nev,evNumber2,nmodules," "," ");
 175      //get pointer to modules array
 176      TObjArray *ITSmodules = ITS->GetModules();
 177      AliITShit *itsHit;
 178
 179      // get the Tree for clusters
 180      ITS->GetTreeC(nev);
 181      TTree *TC=ITS->TreeC();
 182      Int_t nent=TC->GetEntries();
 183      printf("Found %d entries in the tree (must be one per module per event!)\n",nent);
 184      Int_t lay, lad, det;
 185      AliITSgeom *geom = ITS->GetITSgeom();
 186
 187    AliITSDetType *iDetType=ITS->DetType(0);
 188    AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
 189    printf("SPD dimensions %f %f %f \n",seg0->Dx(),seg0->Dz(),seg0->Dy());
 190    printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx());
 191    printf("SPD pitches %d %d \n",seg0->Dpz(0),seg0->Dpx(0));
 192    Float_t SPDlength = seg0->Dz();
 193    Float_t SPDwidth = seg0->Dx();
 194    Float_t SPDthickness = seg0->Dy();
 195    Float_t xpitch = seg0->Dpx(0);
 196    if(SPDlength > 80000) SPDthickness = 150;
 197    Float_t ylim;
 198    if(SPDthickness < 200) {
 199      ylim = SPDthickness/2 - 4;
 200    }else{
 201      ylim = SPDthickness/2 - 10;
 202    }
 203
 204      for (Int_t idettype=0;idettype<3;idettype++) {
 205
 206        TClonesArray *ITSclusters  = ITS->ClustersAddress(idettype);
 207        //printf ("ITSclusters %p \n",ITSclusters);
 208
 209           if (idettype != 0) continue;
 210
 211           Float_t occup1 = 0;
 212           Float_t occup2 = 0;
 213
 214           // Module loop
 215           for (Int_t mod=0; mod<nent; mod++) {
 216               AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod);
 217               geom->GetModuleId(mod,lay,lad,det);
 218
 219               Int_t nhits = itsModule->GetNhits();
 220               //if(nhits) printf("module nhits %d %d\n",mod,nhits);
 221               if(!nhits) continue;
 222
 223               ITS->ResetClusters();
 224               TC->GetEvent(mod);
 225               Int_t nclust = ITSclusters->GetEntries();
 226               if (!nclust) continue;
 227
 228               // cluster/hit loops
 229               //cout<<"mod,lay,nclust,nhits ="<<mod<<","<<lay<<","<<nclust<<","<<nhits<<endl;
 230         for (Int_t clu=0;clu<nclust;clu++) {
 231                 itsclu   = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu);
 232
 233                 Int_t noverlaps = 0;
 234                 Int_t noverprim = 0;
 235
 236                 Int_t clustersizex = itsclu->NclX();
 237                 Int_t clustersizez = itsclu->NclZ();
 238                 Int_t xstart = itsclu->XStartf();
 239                 Int_t xstop = itsclu->XStopf();
 240                 Float_t fxstart = xstart*xpitch;
 241                 Float_t fxstop = (xstop+1)*xpitch;
 242                 Float_t zstart = itsclu->ZStart();
 243                 Float_t zstop = itsclu->ZStop();
 244                 Int_t zend = itsclu->Zend();
 245                 Int_t ntrover = itsclu->NTracks();
 246                 Float_t clusterx = 0;
 247                 Float_t clusterz = itsclu->Z();
 248
 249                 Int_t tr0, tr1, tr2;
 250                 itsclu->GetTracks(tr0,tr1,tr2);
 251
 252                 for(Int_t ii=0;ii<clustersizex;ii++) {
 253                   clusterx += (xstart+0.5+ii)*xpitch;
 254                 }
 255                 clusterx /= clustersizex;
 256                 clusterz /= clustersizez;
 257
 258                 Float_t clusterQ = itsclu->Q();
 259
 260                 if(lay == 1) occup1 += clusterQ;
 261                 if(lay == 2) occup2 += clusterQ;
 262
 263                 ntuple2_st.lay = lay;
 264                 ntuple2_st.x = clusterx/1000.;
 265                 ntuple2_st.z = clusterz/1000.;
 266                 ntuple2_st.nx = clustersizex;
 267                 ntuple2_st.nz = clustersizez;
 268                 ntuple2_st.qcl = clusterQ;
 269
 270                 ntuple2->Fill();
 271
 272                 Int_t icl = 0;
 273                 Float_t dxprimlast = 10.e+6;
 274                 Float_t dzprimlast = 10.e+6;
 275
 276                 Float_t xhit0 = 1e+6;
 277                 Float_t zhit0 = 1e+6;
 278        for (Int_t hit=0;hit<nhits;hit++) {
 279
 280                   // Find coordinate differences between the hit and cluster positions
 281                   // for the resolution determination.
 282
 283                   itsHit   = (AliITShit*)itsModule->GetHit(hit);
 284
 285                   Int_t hitlayer = itsHit->GetLayer();
 286                   Int_t hitladder= itsHit->GetLadder();
 287                   Int_t hitdet= itsHit->GetDetector();
 288                   Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
 289                   Int_t track = itsHit->GetTrack();
 290                   Int_t dray = 0;
 291                   Int_t hitstat = itsHit->GetTrackStatus();
 292                   Float_t zhit = 10000*itsHit->GetZL();
 293                   Float_t xhit = 10000*itsHit->GetXL();
 294                   Float_t yhit = 10000*itsHit->GetYL();
 295         Int_t parent = itsHit->GetParticle()->GetFirstMother();
 296         Int_t partcode = itsHit->GetParticle()->GetPdgCode();
 297
 298         Int_t hitprim = 0;
 299
 300         if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
 301         // particles
 302
 303         Float_t pxsimL = itsHit->GetPXL();  // the momenta at GEANT points
 304         Float_t pysimL = itsHit->GetPYL();
 305         Float_t pzsimL = itsHit->GetPZL();
 306         Float_t psimL = TMath::Sqrt(pxsimL*pxsimL+pysimL*pysimL+pzsimL*pzsimL);
 307
 308         // Check boundaries
 309         if(zhit  > SPDlength/2) {
 310           //cout<<"!!! z outside ="<<zhit<<endl;
 311          zhit = SPDlength/2 - 10;
 312         }
 313         if(zhit < 0 && zhit < -SPDlength/2) {
 314           //cout<<"!!! z outside ="<<zhit<<endl;
 315          zhit = -SPDlength/2 + 10;
 316         }
 317         if(xhit  > SPDwidth/2) {
 318           //cout<<"!!! x outside ="<<xhit<<endl;
 319          xhit = SPDwidth/2 - 10;
 320         }
 321         if(xhit  < 0 && xhit < -SPDwidth/2) {
 322           //cout<<"!!! x outside ="<<xhit<<endl;
 323          xhit = -SPDwidth/2 + 10;
 324         }
 325
 326                   zhit += SPDlength/2;
 327                   xhit += SPDwidth/2;
 328
 329
 330
 331                   if(hitlayer == 1 && hitstat == 66 && yhit > ylim) {
 332                     //if(hitstat == 66 && hitprim ==1) {
 333                     xhit0 = xhit;
 334                     zhit0 = zhit;
 335                   }
 336
 337
 338                   if(hitlayer == 2 && hitstat == 66 && yhit < -ylim) {
 339                     xhit0 = xhit;
 340                     zhit0 = zhit;
 341                   }
 342
 343                   if(hitstat != 68) continue; // Take only the hit if the last
 344                   // track point went out from
 345                   // the detector.
 346
 347                   if(xhit0 > 1e+5 || zhit0 > 1e+5) continue;
 348
 349                   Float_t xmed = (xhit + xhit0)/2;
 350                   Float_t zmed = (zhit + zhit0)/2;
 351
 352                   Float_t xdif = xmed - clusterx;
 353                   Float_t zdif = zmed - clusterz;
 354                   xhit0 = 1e+6;
 355                   zhit0 = 1e+6;
 356
 357         // Consider the hits inside of cluster region only   b.b.
 358
 359                   if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
 360
 361           // Consider the hits only with the track number equaled to one
 362           // of the cluster
 363                     //if((track == tr0) || (track == tr1) || (track == tr2)) {
 364
 365         icl = 1;
 366
 367                 //        part = (TParticle *)particles.UncheckedAt(track);
 368                 //        Int_t partcode = part->GetPdgCode();
 369                 //              Int_t primery = gAlice->GetPrimary(track);
 370
 371
 372 //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
 373 //                      310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 374
 375         Float_t pmod = itsHit->GetParticle()->P(); // total momentum at the
 376                                                    // vertex
 377         Float_t energy = itsHit->GetParticle()->Energy(); // energy at the
 378                                                    // vertex
 379         Float_t mass = itsHit->GetParticle()->GetMass(); // particle mass
 380
 381         Float_t pz = itsHit->GetParticle()->Pz(); // z momentum componetnt
 382                                                    // at the vertex
 383         Float_t px = itsHit->GetParticle()->Px(); // z momentum componetnt
 384                                                    // at the vertex
 385         Float_t py = itsHit->GetParticle()->Py(); // z momentum componetnt
 386                                                    // at the vertex
 387         Float_t phi = itsHit->GetParticle()->Phi(); // Phi angle at the
 388                                                    // vertex
 389         Float_t theta = itsHit->GetParticle()->Theta(); // Theta angle at the
 390                                                    // vertex
 391         //Float_t eta = itsHit->GetParticle()->Eta(); // Pseudo rapidity at the
 392                                                    // vertex
 393         Float_t y;
 394         if((energy-pz) > 0) {
 395           y = 0.5*TMath::Log((energy+pz)/(energy-pz));
 396         }else{
 397           cout<<" Warning: energy < pz ="<<energy<<","<<pz<<endl;
 398           y = 10;
 399         }
 400         Float_t eta = -TMath::Log(TMath::Tan(theta/2));
 401         pmod *= 1.0e+3;
 402
 403
 404         Float_t pxsim = itsHit->GetPXG();  // the momenta at this GEANT point
 405         Float_t pysim = itsHit->GetPYG();
 406         Float_t pzsim = itsHit->GetPZG();
 407         Float_t psim = TMath::Sqrt(pxsim*pxsim+pysim*pysim+pzsim*pzsim);
 408
 409         //Int_t hitprim = 0;
 410
 411         if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 412                                                  // at p < 6 MeV/c
 413
 414         if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but
 415                                                  // not for delta ray which
 416                                                  // also went out from the
 417                                                  // detector and returned
 418                                                  // again
 419
 420         //if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
 421                                               // particles
 422
 423         if(hitprim > 0) noverprim = noverprim + 1;
 424
 425         if(hitprim > 0) {
 426          dxprimlast = xdif;
 427          dzprimlast = zdif;
 428         }
 429
 430         // fill ntuple
 431
 432          ntuple_st.lay = hitlayer;
 433          ntuple_st.nx = clustersizex;
 434          ntuple_st.nz = clustersizez;
 435          ntuple_st.hitprim = hitprim;
 436          ntuple_st.partcode = partcode;
 437          ntuple_st.ntrover = ntrover;
 438          ntuple_st.dx = xdif;
 439          ntuple_st.dz = zdif;
 440          ntuple_st.pmod = pmod;
 441
 442          ntuple->Fill();
 443
 444         if(hitlayer == 1) {
 445            Y1DenA->Fill(y);
 446            Eta1DenA->Fill(eta);
 447         }
 448         if(hitlayer == 2) {
 449            Y2DenA->Fill(y);
 450            Eta2DenA->Fill(eta);
 451         }
 452
 453
 454
 455       if(hitprim > 0) {   // for primary particles
 456
 457         if(hitlayer == 1) {
 458            Xres1->Fill(xdif);
 459            Zres1->Fill(zdif);
 460            Ptot1->Fill(pmod/1000.);
 461            Pz1->Fill(pz);
 462            Theta1->Fill(theta);
 463            Y1->Fill(y);
 464            Eta1->Fill(eta);
 465            Y1Den->Fill(y);
 466            Eta1Den->Fill(eta);
 467            Phi1->Fill(phi);
 468         }
 469         if(hitlayer == 2) {
 470            Xres2->Fill(xdif);
 471            Zres2->Fill(zdif);
 472            Ptot2->Fill(pmod/1000.);
 473            Pz2->Fill(pz);
 474            Theta2->Fill(theta);
 475            Y2->Fill(y);
 476            Eta2->Fill(eta);
 477            Y2Den->Fill(y);
 478            Eta2Den->Fill(eta);
 479            Phi2->Fill(phi);
 480         }
 481       } // primery particles
 482
 483      } // end of cluster region
 484
 485    } // end of hit loop
 486
 487       if(icl == 1) {
 488
 489         // fill ntuple1
 490
 491         if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
 492                                           // delta rays only
 493
 494          ntuple1_st.lay = lay;
 495          ntuple1_st.lad = lad;
 496          ntuple1_st.det = det;
 497          ntuple1_st.x = clusterx*1000.;
 498          ntuple1_st.z = clusterz*1000.;
 499          ntuple1_st.nx = clustersizex;
 500          ntuple1_st.nz = clustersizez;
 501          ntuple1_st.qcl = clusterQ;
 502          ntuple1_st.ntrover = ntrover;
 503          ntuple1_st.noverlaps = noverlaps;
 504          ntuple1_st.noverprim = noverprim;
 505          ntuple1_st.dx = dxprimlast;
 506          ntuple1_st.dz = dzprimlast;
 507
 508          ntuple1->Fill();
 509
 510      } // icl = 1
 511     } // cluster loop
 512    } // module loop
 513
 514      cout<<" Occupancy for layer-1 ="<<occup1<<endl;
 515      cout<<" Occupancy for layer-2 ="<<occup2<<endl;
 516      // The real occupancy values are:
 517      // (for full ALICE event at the full SPD acceptence)
 518      //   occup1 /= 3932160;
 519      //   occup2 /= 7864320;
 520
 521   } // idettype loop
 522  } // end if ITS
 523 } // event loop
 524
 525
 526    //  Write and Draw Histogramms and ntuples
 527
 528
 529
 530    TFile fhistos("SPD_his_dubna.root","RECREATE");
 531
 532    ntuple->Write();
 533    ntuple1->Write();
 534    ntuple2->Write();
 535
 536    Nxpix1->Write();
 537    Nzpix1->Write();
 538    Nxpix2->Write();
 539    Nzpix2->Write();
 540
 541    Xpix1->Write();
 542    Zpix1->Write();
 543    Xpix2->Write();
 544    Zpix2->Write();
 545
 546    Xres1->Write();
 547    Zres1->Write();
 548    Xres2->Write();
 549    Zres2->Write();
 550
 551    Ptot1->Write();
 552    Pz1->Write();
 553    Theta1->Write();
 554    Y1->Write();
 555    Eta1->Write();
 556    Y1Den->Write();
 557    Eta1Den->Write();
 558    Y1DenA->Write();
 559    Eta1DenA->Write();
 560    Phi1->Write();
 561
 562    Ptot2->Write();
 563    Pz2->Write();
 564    Theta2->Write();
 565    Y2->Write();
 566    Eta2->Write();
 567    Y2Den->Write();
 568    Eta2Den->Write();
 569    Y2DenA->Write();
 570    Eta2DenA->Write();
 571    Phi2->Write();
 572
 573    fhistos.Close();
 574    cout<<"!!! Histogramms and ntuples were written"<<endl;
 575
 576    TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
 577    c1->Divide(2,2);
 578    c1->cd(1);
 579    gPad->SetFillColor(33);
 580          Xres1->SetFillColor(42);
 581          Xres1->Draw();
 582    c1->cd(2);
 583    gPad->SetFillColor(33);
 584          Zres1->SetFillColor(46);
 585          Zres1->Draw();
 586    c1->cd(3);
 587    gPad->SetFillColor(33);
 588          Xres2->SetFillColor(42);
 589          Xres2->Draw();
 590    c1->cd(4);
 591    gPad->SetFillColor(33);
 592          Zres2->SetFillColor(46);
 593          Zres2->Draw();
 594
 595      cout<<"END  test for clusters and hits "<<endl;
 596
 597 //     file->Close();
 598 }
 599
 600
 601