ITS/SSDrecpointTest.C

   1 void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
   2   //void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=999)
   3 {
   4 /////////////////////////////////////////////////////////////////////////
   5 //   This macro is a small example of a ROOT macro
   6 //   illustrating how to read the output of GALICE
   7 //   and fill some histograms.
   8 //
   9 //     Root > .L anal.C   //this loads the macro in memory
  10 //     Root > anal();     //by default process first event
  11 //     Root > anal(2);    //process third event
  12 //Begin_Html
  13 /*
  14 <img src="gif/anal.gif">
  15 */
  16 //End_Html
  17 /////////////////////////////////////////////////////////////////////////
  18
  19 // Dynamically link some shared libs
  20
  21    if (gClassTable->GetID("AliRun") < 0) {
  22       gROOT->LoadMacro("loadlibs.C");
  23       loadlibs();
  24    }
  25
  26 // Connect the Root Galice file containing Geometry, Kine and Hits
  27    TString *str = new TString("galice.root");
  28    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(str->Data());
  29    if (!file) file = new TFile(str->Data(),"UPDATE");
  30
  31 // Get AliRun object from file or create it if not on file
  32    //   if (!gAlice) {
  33      gAlice = (AliRun*)file->Get("gAlice");
  34      if (gAlice) printf("AliRun object found on file\n");
  35      if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  36      //}
  37
  38
  39      // -------------- Create ntuples --------------------
  40
  41      //  ntuple structures:
  42
  43
  44           struct {
  45             Int_t lay;
  46             Int_t nxP;
  47             Int_t nxN;
  48             Int_t hitprim;
  49             Int_t partcode;
  50             Int_t ntrover;
  51             Float_t x;
  52             Float_t z;
  53             Float_t dx;
  54             Float_t dz;
  55             Float_t pmod;
  56           } ntuple_st;
  57
  58           struct {
  59             Int_t lay;
  60             Int_t lad;
  61             Int_t det;
  62             Int_t nxP;
  63             Int_t nxN;
  64             Int_t noverlaps;
  65             Int_t noverprim;
  66             Int_t ntrover;
  67             Float_t qclP;
  68             Float_t qclN;
  69             Float_t qrec;
  70             Float_t dx;
  71             Float_t dz;
  72           } ntuple1_st;
  73
  74           struct {
  75             Int_t nxP;
  76             Int_t nxN;
  77             Float_t x;
  78             Float_t z;
  79           } ntuple2_st;
  80
  81           ntuple = new TTree("ntuple","Demo ntuple");
  82           ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
  83           ntuple->Branch("nxP",&ntuple_st.nxP,"nxP/I");
  84           ntuple->Branch("nxN",&ntuple_st.nxN,"nxN/I");
  85           ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
  86           ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
  87           ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
  88           ntuple->Branch("x",&ntuple_st.x,"x/F");
  89           ntuple->Branch("z",&ntuple_st.z,"z/F");
  90           ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
  91           ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
  92           ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
  93
  94           ntuple1 = new TTree("ntuple1","Demo ntuple1");
  95           ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
  96           ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
  97           ntuple1->Branch("det",&ntuple1_st.det,"det/I");
  98           ntuple1->Branch("nxP",&ntuple1_st.nxP,"nxP/I");
  99           ntuple1->Branch("nxN",&ntuple1_st.nxN,"nxN/I");
 100           ntuple1->Branch("qclP",&ntuple1_st.qclP,"qclP/F");
 101           ntuple1->Branch("qclN",&ntuple1_st.qclN,"qclN/F");
 102           ntuple1->Branch("qrec",&ntuple1_st.qrec,"qrec/F");
 103           ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
 104           ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
 105           ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
 106           ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
 107           ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
 108
 109           ntuple2 = new TTree("ntuple2","Demo ntuple2");
 110           ntuple2->Branch("nxP",&ntuple2_st.nxP,"nxP/I");
 111           ntuple2->Branch("nxN",&ntuple2_st.nxN,"nxN/I");
 112           ntuple2->Branch("x",&ntuple2_st.x,"x/F");
 113           ntuple2->Branch("z",&ntuple2_st.z,"z/F");
 114
 115
 116           // Create Histogramms
 117
 118           TH1F *NxP5 = new TH1F("NxP5","P cluster size for layer 5",20,0.,20.);
 119           TH1F *NxN5 = new TH1F("NxN5","N cluster size for layer 5",20,0.,20.);
 120           TH1F *NxP6 = new TH1F("NxP6","P cluster size for layer 6",20,0.,20.);
 121           TH1F *NxN6 = new TH1F("NxN6","N cluster size for layer 6",20,0.,20.);
 122
 123           TH1F *Xres5 = new TH1F("Xres5","Xrec and Xgen difference (micr) for layers 5",100,-200.,200.);
 124           TH1F *Xres6 = new TH1F("Xres6","Xrec and Xgen difference (micr) for layers 6",100,-200.,200.);
 125           TH1F *Zres5 = new TH1F("Zres5","Zrec and Zgen difference (micr) for layers 5",100,-8000.,8000.);
 126           TH1F *Zres6 = new TH1F("Zres6","Zrec and Zgen difference (micr) for layers 6",100,-8000.,8000.);
 127           TH1F *Path5 = new TH1F("Path5","Path length in Si",100,0.,600.);
 128           TH1F *Path6 = new TH1F("Path6","Path length in Si",100,0.,600.);
 129           TH1F *dEdX = new TH1F("dEdX","dEdX  (KeV)",100,0.,500.);
 130           TH2F *adcPadcN5all = new TH2F("adcPadcN5all","adcP/N correlation for lay5",100,0.,200.,100,0.,200.);
 131           TH2F *adcPadcN6all = new TH2F("adcPadcN6all","adcP/N correlation for lay6",100,0.,200.,100,0.,200.);
 132           TH2F *adcPadcN5cut = new TH2F("adcPadcN5cut","adcP/N correlation for lay5 and cut of P-N signas",100,0.,200.,100,0.,200.);
 133           TH2F *adcPadcN6cut = new TH2F("adcPadcN6cut","adcP/N correlation for lay6 and cut of P-N signals",100,0.,200.,100,0.,200.);
 134
 135
 136    AliITS *ITS  = (AliITS*) gAlice->GetModule("ITS");
 137    if (!ITS) { cout << "no ITS" << endl; return; }
 138
 139    //AliITSgeom *aliitsgeo = ITS->GetITSgeom();
 140    AliITSgeom *geom = ITS->GetITSgeom();
 141
 142    //Int_t cp[8]={0,0,0,0,0,0,0,0};
 143
 144    cout << "SSD" << endl;
 145
 146    AliITSDetType *iDetType=ITS->DetType(2);
 147    AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
 148    AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
 149    //res2->SetSigmaSpread(3.,2.);
 150    AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
 151    ITS->SetSimulationModel(2,sim2);
 152
 153    TClonesArray *dig2  = ITS->DigitsAddress(2);
 154    TClonesArray *recp2  = ITS->ClustersAddress(2);
 155    //   AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2,recp2);
 156    AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
 157    ITS->SetReconstructionModel(2,rec2);
 158    // test
 159    printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz());
 160    printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
 161
 162
 163 //
 164 //   Loop over events
 165 //
 166
 167
 168    Int_t Nh=0;
 169    Int_t Nh1=0;
 170    for (int nev=0; nev<= evNumber2; nev++) {
 171      Int_t nparticles = 0;
 172      nparticles = gAlice->GetEvent(nev);
 173      cout << "nev         " << nev <<endl;
 174      cout << "nparticles  " << nparticles <<endl;
 175      if (nev < evNumber1) continue;
 176      if (nparticles <= 0) return;
 177
 178      AliITShit *itsHit;
 179      AliITSRecPoint *itsPnt = 0;
 180      AliITSRawClusterSSD *itsClu = 0;
 181
 182      // Get Hit, Cluster & Recpoints Tree Pointers
 183
 184      TTree *TH = gAlice->TreeH();
 185      Int_t nenthit=TH->GetEntries();
 186      printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
 187
 188      ITS->GetTreeC(nev);
 189      TTree *TC=ITS->TreeC();
 190      Int_t nentclu=TC->GetEntries();
 191      printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);
 192
 193      TTree *TR = gAlice->TreeR();
 194      Int_t nentrec=TR->GetEntries();
 195      printf("Found %d entries in the RecPoints tree\n",nentrec);
 196
 197      // Get Pointers to Clusters & Recpoints TClonesArrays
 198
 199      TClonesArray *ITSclu  = ITS->ClustersAddress(2);
 200      printf ("ITSclu %p \n",ITSclu);
 201      TClonesArray *ITSrec  = ITS->RecPoints();
 202      printf ("ITSrec %p \n",ITSrec);
 203
 204      // check recpoints
 205
 206      //Int_t nbytes = 0;
 207      Int_t totpoints = 0;
 208      Int_t totclust = 0;
 209
 210      // check hits
 211
 212      Int_t nmodules=0;
 213      Int_t mod;
 214
 215      ITS->InitModules(-1,nmodules);
 216      ITS->FillModules(nev,0,nmodules,"","");
 217
 218      TObjArray *fITSmodules = ITS->GetModules();
 219
 220      Int_t first0 = geom->GetStartDet(0);  // SPD
 221      Int_t last0 = geom->GetLastDet(0);    // SPD
 222      Int_t first1 = geom->GetStartDet(1);  // SDD
 223      Int_t last1 = geom->GetLastDet(1);    // SDD
 224      Int_t first2 = geom->GetStartDet(2);  // SSD
 225      Int_t last2 = geom->GetLastDet(2);    // SSD
 226
 227      //  For the SPD: first0 = 0, last0 = 239     (240 modules);
 228      //  for the SDD: first1 = 240, last1 = 499   (260 modules);
 229      //  for the SSD: first2 = 500, last2 = 2269  (1770 modules).
 230
 231      printf("det type %d first0, last0 %d %d \n",0,first0,last0);
 232      printf("det type %d first1, last1 %d %d \n",1,first1,last1);
 233      printf("det type %d first2, last2 %d %d \n",2,first2,last2);
 234
 235      // module loop for the SSD
 236      //for (mod=first2; mod<last2+1; mod++) {  // for the "ALL" option
 237      for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option
 238
 239        TTree *TR = gAlice->TreeR();
 240        Int_t nentrec=TR->GetEntries();
 241        //printf("Found %d entries in the RecPoints tree\n",nentrec);
 242
 243               //cout << "CLUSTERS: reset" << endl;
 244        ITS->ResetClusters();
 245        //cout << "CLUSTERS: get" << endl;
 246        TC->GetEvent(mod);
 247        //cout << "RECPOINTS: reset" << endl;
 248        ITS->ResetRecPoints();
 249        //cout << "RECPOINTS: get" << endl;
 250        //TR->GetEvent(mod+1);   // for the V3.04 AliRoot
 251        TR->GetEvent(mod);       // for the V3.05 AliRoot
 252
 253        Int_t nrecp = ITSrec->GetEntries();
 254        totpoints += nrecp;
 255        if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
 256        //if (!nrecp) continue;
 257        Int_t nclusters = ITSclu->GetEntries();
 258        totclust += nclusters;
 259        //if (nclusters) printf("Found %d clusters for module %d\n",nrecc,mod);
 260
 261        AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod+first2);
 262        // for the "SSD" option
 263
 264        //AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
 265        // for the "ALL" option
 266
 267        //       printf("Mod: %X\n",Mod);
 268        Int_t nhits = Mod->GetNhits();
 269        Float_t epart = 0;
 270        //cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
 271
 272        // ---------------- cluster/hit analysis ---------------------
 273
 274
 275      Float_t pathInSSD = 300.;
 276
 277        // ---- Recpoint loop
 278        for (Int_t pnt=0;pnt<nrecp;pnt++) {
 279          itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
 280          if(!itsPnt) continue;
 281          itsClu  = (AliITSRawClusterSSD*)ITSclu->At(pnt);
 282          if(!itsClu) continue;
 283
 284          Int_t nxP = itsClu->fMultiplicity;
 285          Int_t nxN = itsClu->fMultiplicityN;
 286          Int_t ntrover = itsClu->fNtracks;
 287          Float_t qclP = itsClu->fSignalP;     // in ADC
 288          Float_t qclN = itsClu->fSignalN;     // in ADC
 289          //Float_t dq = qclP - qclN;
 290          Float_t qcut = itsClu->fQErr;        // abs(dq)/signal,
 291                                               // where signal is
 292                                               // max of qclP,qclN
 293          Float_t xrec = 10000*itsPnt->GetX();
 294          Float_t zrec = 10000*itsPnt->GetZ();
 295          Float_t qrec = itsPnt->GetQ();      // in ADC, maximum from fSignalP/N
 296          //Float_t dedx = itsPnt->GetdEdX();   // in KeV (ADC * 2.16)
 297          Float_t dedx = itsPnt->fdEdX;   // in KeV (ADC * 2.16)
 298          Int_t ii = 0;
 299          Int_t tr1 = itsPnt->GetLabel(ii);
 300          Int_t ii = 1;
 301          Int_t tr2 = itsPnt->GetLabel(ii);
 302          Int_t ii = 2;
 303          Int_t tr3 = itsPnt->GetLabel(ii);
 304          //cout<<"pnt,ntrover,tr1,tr2,tr3 ="<<pnt<<","<<ntrover<<","<<tr1<<","<<tr2<<","<<tr3<<endl;
 305          // fill ntuple2
 306              ntuple2_st.nxP = nxP;
 307              ntuple2_st.nxN = nxN;
 308              ntuple2_st.x = xrec/1000;
 309              ntuple2_st.z = zrec/1000;
 310
 311              if(qcut < 0.18) ntuple2->Fill();
 312
 313
 314           Int_t noverlaps = 0;
 315           Int_t noverprim = 0;
 316           Int_t flaghit = 0;
 317           Float_t xhit0 = 1e+7;
 318           Float_t yhit0 = 1e+7;
 319           Float_t zhit0 = 1e+7;
 320
 321        // Hit loop
 322         for (Int_t hit=0;hit<nhits;hit++) {
 323
 324          itsHit   = (AliITShit*)Mod->GetHit(hit);
 325
 326          Int_t flagtrack = 0;
 327          Int_t hitlayer = itsHit->GetLayer();
 328          Int_t hitladder= itsHit->GetLadder();
 329          Int_t hitdet= itsHit->GetDetector();
 330
 331          Int_t track = itsHit->GetTrack();
 332          Int_t dray = 0;
 333          Int_t hitstat = itsHit->GetTrackStatus();
 334
 335           Float_t zhit = 10000*itsHit->GetZL();
 336           Float_t xhit = 10000*itsHit->GetXL();
 337           Float_t yhit = 10000*itsHit->GetYL();
 338           Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
 339
 340            Int_t parent = itsHit->GetParticle()->GetFirstMother();
 341            Int_t partcode = itsHit->GetParticle()->GetPdgCode();
 342
 343    //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
 344    //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 345
 346            //cout<<"hit,hitstat,track,partcode,parent ="<<hit<<","<<hitstat<<","<<track<<","<<partcode<<","<<parent<<endl;
 347            Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
 348                                                       // vertex
 349            pmod *= 1.0e+3;
 350
 351           if(hitstat == 66 && yhit < -146.) {  // entering hit
 352             xhit0 = xhit;
 353             yhit0 = yhit;
 354             zhit0 = zhit;
 355           }
 356
 357           if(hitstat == 66) continue; // Take the not entering hits only
 358
 359           if(xhit0 > 9e+6 || zhit0 > 9e+6 || yhit0 > 9e+6) {
 360             //cout<<"default xhit0,zhit0,yhit0 ="<<xhit0<<","<<zhit0<<","<<yhit0<<endl;
 361             continue;
 362           }
 363
 364
 365
 366           // Consider the hits only with the track number equaled to one
 367           // of the recpoint
 368           if((track == tr1) || (track == tr2) || (track == tr3)) flagtrack = 1;
 369           //cout<<"!!Ok: track,tr1,tr2,tr3 ="<<track<<","<<tr1<<","<<tr2<<","<<tr3<<endl;
 370
 371          if(flagtrack == 1) {     // the hit corresponds to the recpoint
 372
 373            flaghit = 1;
 374
 375            //Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
 376            //Float_t py = itsHit->GetPYL();
 377            //Float_t pz = itsHit->GetPZL();
 378
 379          Int_t hitprim = 0;
 380
 381          if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 382                                                   // at p < 6 MeV/c
 383
 384          if((hitstat == 68 || hitstat == 33) && dray == 0)  noverlaps=noverlaps + 1;
 385                                                   // overlapps for all hits but
 386                                                   // not for delta ray which
 387                                                   // also went out from the
 388                                                   // detector and returned
 389                                                   // again
 390
 391
 392           // x,z resolution colculation
 393           if((hitstat == 68 || hitsat == 33) && dray == 0) {
 394              Float_t xmed = (xhit + xhit0)/2;
 395              Float_t zmed = (zhit + zhit0)/2;
 396              Float_t xdif = xmed - xrec;
 397              Float_t zdif = zmed - zrec;
 398
 399             if(parent < 0)  {
 400               hitprim = 1; // hitprim=1 for the primery particles
 401               noverprim += 1;
 402             }
 403              pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
 404
 405              //cout<<"lay,pnt,hit,xmed,xrec,xdif,zmed,zrec,zdif ="<<hitlayer<<","<<pnt<<","<<hit<<","<<xmed<<","<<xrec<<","<<xdif<<","<<zmed<<","<<zrec<<","<<zdif<<endl;
 406
 407          // fill ntuple
 408              ntuple_st.lay = hitlayer;
 409              ntuple_st.nxP = nxP;
 410              ntuple_st.nxN = nxN;
 411              ntuple_st.hitprim = hitprim;
 412              ntuple_st.partcode = partcode;
 413              ntuple_st.ntrover = ntrover;
 414              ntuple_st.x = xrec/1000;
 415              ntuple_st.z = zrec/1000;
 416              ntuple_st.dx = xdif;
 417              ntuple_st.dz = zdif;
 418              ntuple_st.pmod = pmod;
 419
 420              //if(qcut < 0.18) ntuple->Fill();
 421              ntuple->Fill();
 422
 423              //if(hitlayer == 5 && qcut < 0.18) {
 424             if(hitlayer == 5 ) {
 425              Xres5->Fill(xdif);
 426              Zres5->Fill(zdif);
 427              Path5->Fill(pathInSSD);
 428             }
 429             //if(hitlayer == 6 && qcut < 0.18) {
 430             if(hitlayer == 6) {
 431              Xres6->Fill(xdif);
 432              Zres6->Fill(zdif);
 433              Path6->Fill(pathInSSD);
 434             }
 435           } // hitstat 68/33
 436          } else {       // non correspondent hit
 437           xhit0 = 1e+7;
 438           zhit0 = 1e+7;
 439          } // end of hit-recpoint correspondence
 440         } // hit loop
 441
 442         if(flaghit == 1) {
 443
 444           if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
 445           // delta rays only
 446
 447           // fill ntuple1
 448           ntuple1_st.lay = hitlayer;
 449           ntuple1_st.lad = hitladder;
 450           ntuple1_st.det = hitdet;
 451           ntuple1_st.nxP = nxP;
 452           ntuple1_st.nxN = nxN;
 453           ntuple1_st.qclP = qclP*300/pathInSSD;
 454           ntuple1_st.qclN = qclN*300/pathInSSD;
 455           ntuple1_st.qrec = qrec*300/pathInSSD;
 456           ntuple1_st.dx = xdif;
 457           ntuple1_st.dz = zdif;
 458           noverlaps -= 1;
 459           noverprim -= 1;
 460           ntuple1_st.noverlaps = noverlaps;
 461           ntuple1_st.noverprim = noverprim;
 462           ntuple1_st.ntrover = ntrover;
 463
 464           //if(qcut < 0.18) ntuple1->Fill();
 465           ntuple1->Fill();
 466
 467           Float_t de = dedx*300./pathInSSD;
 468           dEdX->Fill(de);
 469             if(hitlayer == 5 ) {
 470              adcPadcN5all->Fill(qclP,qclN);
 471             }
 472             if(hitlayer == 6 ) {
 473              adcPadcN6all->Fill(qclP,qclN);
 474             }
 475             if(hitlayer == 5 && qcut < 0.18) {
 476              adcPadcN5cut->Fill(qclP,qclN);
 477              NxP5->Fill(nxP);
 478              NxN5->Fill(nxN);
 479             }
 480             if(hitlayer == 6 && qcut < 0.18) {
 481              adcPadcN6cut->Fill(qclP,qclN);
 482              NxP6->Fill(nxP);
 483              NxN6->Fill(nxN);
 484             }
 485         } // flaghit = 1
 486        } //b.b. recpoint loop
 487      } //b.b. module loop
 488    } //b.b. evnt loop
 489
 490    TFile fhistos("SSD_his.root","RECREATE");
 491
 492    ntuple->Write();
 493    ntuple1->Write();
 494    ntuple2->Write();
 495    NxP5->Write();
 496    NxN5->Write();
 497    NxP6->Write();
 498    NxN6->Write();
 499    Xres5->Write();
 500    Zres5->Write();
 501    Xres6->Write();
 502    Zres6->Write();
 503    Path5->Write();
 504    Path6->Write();
 505    adcPadcN5all->Write();
 506    adcPadcN6all->Write();
 507    adcPadcN5cut->Write();
 508    adcPadcN6cut->Write();
 509    dEdX->Write();
 510
 511    fhistos.Close();
 512
 513    cout<<"!!! Histogramms and ntuples were written"<<endl;
 514
 515    TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
 516    c1->Divide(2,2);
 517    c1->cd(1);
 518    gPad->SetFillColor(33);
 519          Xres5->SetFillColor(42);
 520          Xres5->Draw();
 521    c1->cd(2);
 522    gPad->SetFillColor(33);
 523          Zres5->SetFillColor(46);
 524          Zres5->Draw();
 525    c1->cd(3);
 526    gPad->SetFillColor(33);
 527          Xres6->SetFillColor(42);
 528          Xres6->Draw();
 529    c1->cd(4);
 530    gPad->SetFillColor(33);
 531          Zres6->SetFillColor(46);
 532          Zres6->Draw();
 533
 534    cout<<"END  test for clusters and hits "<<endl;
 535
 536 }
 537
 538
 539
 540
 541
 542