RICH/RICHpadtest.C

   1 void RICHpadtest (Int_t diaglevel,Int_t evNumber1=0,Int_t evNumber2=0)
   2 {
   3
   4 // Diaglevel
   5 // 1-> Single Ring Hits
   6 // 2-> Single Ring Spectra
   7 // 3-> Single Ring Statistics
   8 // 4-> Single Ring Reconstruction
   9 // 5-> Full Event Hits
  10
  11 /////////////////////////////////////////////////////////////////////////
  12 //   This macro is a small example of a ROOT macro
  13 //   illustrating how to read the output of GALICE
  14 //   and do some analysis.
  15 //
  16 /////////////////////////////////////////////////////////////////////////
  17
  18
  19     Int_t NpadX = 162;                 // number of pads on X
  20     Int_t NpadY = 162;                 // number of pads on Y
  21
  22     Int_t Pad[162][162];
  23     for (Int_t i=0;i<NpadX;i++) {
  24         for (Int_t j=0;j<NpadY;j++) {
  25             Pad[i][j]=0;
  26         }
  27     }
  28     gClassTable->GetID("AliRun");
  29
  30
  31 // Dynamically link some shared libs
  32
  33     if (gClassTable->GetID("AliRun") < 0) {
  34       gROOT->LoadMacro("loadlibs.C");
  35       loadlibs();
  36     }
  37     else {
  38       //delete gAlice;
  39       gAlice = 0;
  40     }
  41
  42     gAlice=0;
  43
  44 // Connect the Root Galice file containing Geometry, Kine and Hits
  45
  46     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  47     if (!file) file = new TFile("galice.root","UPDATE");
  48
  49 // Get AliRun object from file or create it if not on file
  50
  51     if (!gAlice) {
  52       gAlice = (AliRun*)file->Get("gAlice");
  53       if (gAlice) printf("AliRun object found on file\n");
  54       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  55     }
  56     else {
  57       delete gAlice;
  58       gAlice = (AliRun*)file->Get("gAlice");
  59       if (gAlice) printf("AliRun object found on file\n");
  60       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  61     }
  62
  63 //  Create some histograms
  64
  65    Int_t xmin= -NpadX/2;
  66    Int_t xmax=  NpadX/2;
  67    Int_t ymin= -NpadY/2;
  68    Int_t ymax=  NpadY/2;
  69
  70    TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-30,30,150,-30,30);
  71
  72    if (diaglevel == 1)
  73
  74      {
  75        printf("Single Ring Hits\n");
  76        TH2F *feedback = new TH2F("feedback","Feedback hit distribution",150,-30,30,150,-30,30);
  77        TH2F *mip = new TH2F("mip","Mip hit distribution",30,-3,3,30,-3,3);
  78        TH2F *cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-30,30,150,-30,30);
  79        TH2F *h = new TH2F("h","Detector hit distribution",150,-30,30,150,-30,30);
  80        TH1F *hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-30,30);
  81        TH1F *hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-30,30);
  82      }
  83    else
  84      {
  85        printf("Full Event Hits\n");
  86
  87        TH2F *feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
  88        TH2F *mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
  89        TH2F *cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
  90        TH2F *h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300);
  91        TH1F *hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
  92        TH1F *hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
  93
  94      }
  95
  96    TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
  97    TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
  98    TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
  99    TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
 100    TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
 101    TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
 102    TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
 103
 104    TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
 105    TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",200,.5,1);
 106    TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
 107    TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
 108    TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
 109    TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
 110    TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
 111    TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
 112    TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
 113    TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
 114    TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
 115    TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
 116    TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
 117    TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
 118    TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
 119    TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
 120    TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
 121    TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
 122    TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
 123    TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
 124    TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
 125    TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
 126    TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
 127    TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence",100,0,15);
 128    TH1F *Omega = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.5,1);
 129    TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",200,0,15);
 130    TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",200,-180,180);
 131    TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",200,.5,1);
 132    TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
 133
 134
 135 //   Start loop over events
 136
 137    Int_t Nh=0;
 138    Int_t pads=0;
 139    Int_t Nh1=0;
 140    Int_t mothers[80000];
 141    Int_t mothers2[80000];
 142    Float_t mom[3];
 143    Int_t nraw=0;
 144    Int_t phot=0;
 145    Int_t feed=0;
 146    Int_t padmip=0;
 147    for (Int_t i=0;i<100;i++) mothers[i]=0;
 148    for (int nev=0; nev<= evNumber2; nev++) {
 149        Int_t nparticles = gAlice->GetEvent(nev);
 150
 151
 152        //cout<<"nev  "<<nev<<endl;
 153        printf ("\nProcessing event: %d\n",nev);
 154        //cout<<"nparticles  "<<nparticles<<endl;
 155        printf ("Particles       : %d\n",nparticles);
 156        if (nev < evNumber1) continue;
 157        if (nparticles <= 0) return;
 158
 159 // Get pointers to RICH detector and Hits containers
 160
 161        AliRICH *RICH  = (AliRICH*)gAlice->GetDetector("RICH");
 162        Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
 163        gAlice->TreeR()->GetEvent(nent-1);
 164        TClonesArray *Rawclusters = RICH->RawClustAddress(2);    //  Raw clusters branch
 165        //printf ("Rawclusters:%p",Rawclusters);
 166        Int_t nrawclusters = Rawclusters->GetEntriesFast();
 167        //printf (" nrawclusters:%d\n",nrawclusters);
 168        gAlice->TreeR()->GetEvent(nent-1);
 169        TClonesArray *RecHits = RICH->RecHitsAddress(2);
 170        Int_t nrechits = RecHits->GetEntriesFast();
 171        //printf (" nrechits:%d\n",nrechits);
 172        TTree *TH = gAlice->TreeH();
 173        Int_t ntracks = TH->GetEntries();
 174
 175
 176
 177        Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
 178        gAlice->TreeD()->GetEvent(nent-1);
 179
 180 // Start loop on tracks in the hits containers
 181        Int_t Nc=0;
 182        for (Int_t track=0; track<ntracks;track++) {
 183            printf ("Processing Track: %d\n",track);
 184            gAlice->ResetHits();
 185            Int_t nbytes += TH->GetEvent(track);
 186            if (RICH)  {
 187                //RICH->ResetRawClusters();
 188                TClonesArray *PadHits = RICH->PadHits();      // Cluster branch
 189                TClonesArray *Hits = RICH->Hits();            // Hits branch
 190                TClonesArray *Cerenkovs = RICH->Cerenkovs();  // Cerenkovs branch
 191            }
 192            //see hits distribution
 193
 194
 195            Int_t nhits = Hits->GetEntriesFast();
 196            if (nhits) Nh+=nhits;
 197            //printf("nhits %d\n",nhits);
 198            for (Int_t hit=0;hit<nhits;hit++) {
 199               mHit = (AliRICHHit*) Hits->UncheckedAt(hit);
 200               Int_t nch  = mHit->fChamber;              // chamber number
 201               Float_t x  = mHit->fX;                    // x-pos of hit
 202               Float_t y  = mHit->fZ;                    // y-pos
 203               Float_t phi = mHit->fPhi;                 //Phi angle of incidence
 204               Float_t theta = mHit->fTheta;             //Theta angle of incidence
 205               Int_t index = mHit->fTrack;
 206               Int_t particle = mHit->fParticle;
 207               Int_t freon = mHit->fLoss;
 208
 209              hitsX->Fill(x,(float) 1);
 210              hitsY->Fill(y,(float) 1);
 211
 212               //printf("Particle:%d\n",particle);
 213
 214               TParticle *current = (TParticle*)(*gAlice->Particles())[index];
 215               //printf("Particle type: %d\n",current->GetPdgCode());
 216
 217               hitsTheta->Fill(theta,(float) 1);
 218               if (RICH->GetDebugLevel() == -1)
 219                   printf("Theta:%f, Phi:%f\n",theta,phi);
 220
 221               //printf("Debug Level:%d\n",RICH->GetDebugLevel());
 222
 223               if (TMath::Abs(particle) < 50000000)
 224                 {
 225                   mip->Fill(x,y,(float) 1);
 226                   if (current->Energy() - current->GetCalcMass()>1 && freon==1)
 227                     {
 228                       hitsPhi->Fill(phi,(float) 1);
 229                       //hitsTheta->Fill(theta,(float) 1);
 230                       //printf("Theta:%f, Phi:%f\n",theta,phi);
 231                     }
 232                 }
 233
 234               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
 235                 {
 236                   pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
 237                 }
 238               if (TMath::Abs(particle)==2212)
 239                 {
 240                   protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
 241                 }
 242               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
 243                   || TMath::Abs(particle)==311)
 244                 {
 245                   kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
 246                 }
 247               if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
 248                 {
 249                   if (current->Energy() - current->GetCalcMass()>1)
 250                     chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
 251                 }
 252               //printf("Hits:%d\n",hit);
 253               //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
 254               // Fill the histograms
 255               Nh1+=nhits;
 256               h->Fill(x,y,(float) 1);
 257                   //}
 258               //}
 259           }
 260
 261            Int_t ncerenkovs = Cerenkovs->GetEntriesFast();
 262
 263            if (ncerenkovs) {
 264                for (Int_t hit=0;hit<ncerenkovs;hit++) {
 265                    cHit = (AliRICHCerenkov*) Cerenkovs->UncheckedAt(hit);
 266                    Int_t nchamber = cHit->fChamber;     // chamber number
 267                    Int_t index =    cHit->fTrack;
 268                    Int_t pindex =   cHit->fIndex;
 269                    Int_t cx  =      cHit->fX;                // x-position
 270                    Int_t cy  =      cHit->fZ;                // y-position
 271                    Int_t cmother =  cHit->fCMother;      // Index of mother particle
 272                    Int_t closs =    cHit->fLoss;           // How did the particle get lost?
 273                   //printf ("Cerenkov hit, X:%d, Y:%d\n",cx,cy);
 274
 275
 276
 277                    TParticle *current = (TParticle*)(*gAlice->Particles())[index];
 278                    Float_t energyckov = current->Energy();
 279
 280                    if (current->GetPdgCode() == 50000051)
 281                    {
 282                        if (closs==4)
 283                        {
 284                            feedback->Fill(cx,cy,(float) 1);
 285                            feed++;
 286                        }
 287                    }
 288                    if (current->GetPdgCode() == 50000050)
 289                    {
 290
 291                      totalphotonsevent->Fill(ncerenkovs,(float) 1);
 292                      phspectra2->Fill(energyckov*1e9,(float) 1);
 293
 294                      if (closs==4)
 295                        {
 296                          cerenkov->Fill(cx,cy,(float) 1);
 297
 298
 299
 300                          TParticle *MIP = (TParticle*)(*gAlice->Particles())[cmother];
 301                          mipHit = (AliRICHHit*) Hits->UncheckedAt(0);
 302                          mom[0] = current->Px();
 303                          mom[1] = current->Py();
 304                          mom[2] = current->Pz();
 305                          /*mom[0] = cHit->fMomX;
 306                            mom[1] = cHit->fMomZ;
 307                            mom[2] = cHit->fMomY;*/
 308                          Float_t energymip = MIP->Energy();
 309                          Float_t Mip_px = mipHit->fMomX;
 310                          Float_t Mip_py = mipHit->fMomY;
 311                          Float_t Mip_pz = mipHit->fMomZ;
 312                          /*Float_t Mip_px = MIP->Px();
 313                            Float_t Mip_py = MIP->Py();
 314                            Float_t Mip_pz = MIP->Pz();*/
 315
 316
 317
 318                          Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
 319                          Float_t rt = TMath::Sqrt(r);
 320                          Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;
 321                          Float_t Mip_rt = TMath::Sqrt(Mip_r);
 322                          Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt);
 323                          Float_t cherenkov = TMath::ACos(coscerenkov);
 324                          ckovangle->Fill(cherenkov,(float) 1);                           //Cerenkov angle calculus
 325                          //printf("Cherenkov: %f\n",cherenkov);
 326                          Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
 327                          hckphi->Fill(ckphi,(float) 1);
 328
 329
 330                          Float_t mix = MIP->Vx();
 331                          Float_t miy = MIP->Vy();
 332                          Float_t mx = mipHit->fX;
 333                          Float_t my = mipHit->fZ;
 334                          //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
 335                          Float_t dx = cx - mx;
 336                          Float_t dy = cy - my;
 337                          //printf("Dx:%f, Dy:%f\n",dx,dy);
 338                          Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
 339                          //printf("Final radius:%f\n",final_radius);
 340                          radius->Fill(final_radius,(float) 1);
 341
 342                          phspectra1->Fill(energyckov*1e9,(float) 1);
 343                          phot++;
 344                        }
 345                      for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
 346                        if (cmother == nmothers){
 347                          if (closs == 4)
 348                            mothers2[cmother]++;
 349                          mothers[cmother]++;
 350                        }
 351                      }
 352                    }
 353                }
 354            }
 355
 356            if (nrawclusters) {
 357                for (Int_t hit=0;hit<nrawclusters;hit++) {
 358                    rcHit = (AliRICHRawCluster*) Rawclusters->UncheckedAt(hit);
 359                    //Int_t nchamber = rcHit->fChamber;     // chamber number
 360                    //Int_t nhit = cHit->fHitNumber;        // hit number
 361                    Int_t qtot = rcHit->fQ;                 // charge
 362                    Int_t fx  =  rcHit->fX;                 // x-position
 363                    Int_t fy  =  rcHit->fY;                 // y-position
 364                    Int_t type = rcHit->fCtype;             // cluster type ?
 365                    Int_t mult = rcHit->fMultiplicity;      // How many pads form the cluster
 366                    pads += mult;
 367                    if (qtot > 0) {
 368                        if (fx>-4 && fx<4 && fy>-4 && fy<4) {
 369                            padmip+=mult;
 370                        } else {
 371                            padnumber->Fill(mult,(float) 1);
 372                            nraw++;
 373                            if (mult<4) Clcharge->Fill(qtot,(float) 1);
 374                        }
 375                    }
 376                }
 377            }
 378
 379            if(nrechits)
 380              {
 381                for (Int_t hit=0;hit<nrechits;hit++) {
 382                  recHit = (AliRICHRecHit*) RecHits->UncheckedAt(hit);
 383                  Float_t r_omega = recHit->fOmega;                  // Cerenkov angle
 384                  Float_t r_theta = recHit->fTheta;                  // Theta angle of incidence
 385                  Float_t r_phi   = recHit->fPhi;                    // Phi angle if incidence
 386                  Float_t *cer_pho = recHit->fCerPerPhoton;        // Cerenkov angle per photon
 387                  Int_t *padsx = recHit->fPadsUsedX;           // Pads Used fo reconstruction (x)
 388                  Int_t *padsy = recHit->fPadsUsedY;           // Pads Used fo reconstruction (y)
 389                  Int_t goodPhotons = recHit->fGoodPhotons;    // Number of pads used for reconstruction
 390
 391                  Omega->Fill(r_omega,(float) 1);
 392                  Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
 393                  Phi->Fill(r_phi*180/TMath::Pi(),(float) 1);
 394
 395                  for (Int_t i=0; i<goodPhotons; i++)
 396                    {
 397                      PhotonCer->Fill(cer_pho[i],(float) 1);
 398                      PadsUsed->Fill(padsx[i],padsy[i],1);
 399                      //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]);
 400                    }
 401
 402                  //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
 403                }
 404              }
 405        }
 406
 407        for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
 408            totalphotonstrack->Fill(mothers[nmothers],(float) 1);
 409            mother->Fill(mothers2[nmothers],(float) 1);
 410            //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]);
 411        }
 412
 413        clusev->Fill(nraw,(float) 1);
 414        photev->Fill(phot,(float) 1);
 415        feedev->Fill(feed,(float) 1);
 416        padsmip->Fill(padmip,(float) 1);
 417        padscl->Fill(pads,(float) 1);
 418        //printf("Photons:%d\n",phot);
 419        phot = 0;
 420        feed = 0;
 421        pads = 0;
 422        nraw=0;
 423        padmip=0;
 424
 425        TClonesArray *Digits = RICH->DigitsAddress(2);    //  Raw clusters branch
 426        Int_t ndigits = Digits->GetEntriesFast();
 427        //printf("Digits:%d\n",ndigits);
 428        padsev->Fill(ndigits,(float) 1);
 429
 430        for (Int_t ich=0;ich<7;ich++)
 431          {
 432            TClonesArray *Digits = RICH->DigitsAddress(ich);    //  Raw clusters branch
 433            Int_t ndigits = Digits->GetEntriesFast();
 434            //printf("Digits:%d\n",ndigits);
 435            padsev->Fill(ndigits,(float) 1);
 436            if (ndigits) {
 437              for (Int_t hit=0;hit<ndigits;hit++) {
 438                dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
 439                //Int_t nchamber = padHit->fChamber;     // chamber number
 440                //Int_t nhit = dHit->fHitNumber;          // hit number
 441                Int_t qtot = dHit->fSignal;                // charge
 442                Int_t ipx  = dHit->fPadX;               // pad number on X
 443                Int_t ipy  = dHit->fPadY;               // pad number on Y
 444                //Int_t iqpad  = dHit->fQpad;           // charge per pad
 445                //Int_t rpad  = dHit->fRSec;            // R-position of pad
 446                //printf ("Pad hit, PadX:%d, PadY:%d\n",ipx,ipy);
 447                if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot);
 448                if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
 449                if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
 450                if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
 451                if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
 452                if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
 453                if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
 454                if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
 455              }
 456            }
 457          }
 458    }
 459
 460
 461    //Create canvases, set the view range, show histograms
 462
 463    switch(diaglevel)
 464      {
 465      case 1:
 466
 467        TCanvas *c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350);
 468        hc0->SetXTitle("ix (npads)");
 469        hc0->Draw("box");
 470
 471 //
 472        TCanvas *c4 = new TCanvas("c4","Hits per type",100,100,600,700);
 473        c4->Divide(2,2);
 474
 475        c4->cd(1);
 476        feedback->SetXTitle("x (cm)");
 477        feedback->SetYTitle("y (cm)");
 478        feedback->Draw();
 479
 480        c4->cd(2);
 481        //mip->SetFillColor(42);
 482        mip->SetXTitle("x (cm)");
 483        mip->SetYTitle("y (cm)");
 484        mip->Draw();
 485
 486        c4->cd(3);
 487        //cerenkov->SetFillColor(42);
 488        cerenkov->SetXTitle("x (cm)");
 489        cerenkov->SetYTitle("y (cm)");
 490        cerenkov->Draw();
 491
 492        c4->cd(4);
 493        //h->SetFillColor(42);
 494        h->SetXTitle("x (cm)");
 495        h->SetYTitle("y (cm)");
 496        h->Draw();
 497
 498        TCanvas *c10 = new TCanvas("c10","Hits distribution",150,150,600,350);
 499        c10->Divide(2,1);
 500
 501        c10->cd(1);
 502        hitsX->SetFillColor(42);
 503        hitsX->SetXTitle("(cm)");
 504        hitsX->Draw();
 505
 506        c10->cd(2);
 507        hitsY->SetFillColor(42);
 508        hitsY->SetXTitle("(cm)");
 509        hitsY->Draw();
 510
 511
 512        break;
 513 //
 514      case 2:
 515
 516        TCanvas *c6 = new TCanvas("c6","Photon Spectra",50,50,600,350);
 517        c6->Divide(2,1);
 518
 519        c6->cd(1);
 520        phspectra2->SetFillColor(42);
 521        phspectra2->SetXTitle("energy (eV)");
 522        phspectra2->Draw();
 523        c6->cd(2);
 524        phspectra1->SetFillColor(42);
 525        phspectra1->SetXTitle("energy (eV)");
 526        phspectra1->Draw();
 527
 528        TCanvas *c9 = new TCanvas("c9","Particles Spectra",100,100,600,700);
 529        c9->Divide(2,2);
 530
 531        c9->cd(1);
 532        pionspectra->SetFillColor(42);
 533        pionspectra->SetXTitle("(GeV)");
 534        pionspectra->Draw();
 535
 536        c9->cd(2);
 537        protonspectra->SetFillColor(42);
 538        protonspectra->SetXTitle("(GeV)");
 539        protonspectra->Draw();
 540
 541        c9->cd(3);
 542        kaonspectra->SetFillColor(42);
 543        kaonspectra->SetXTitle("(GeV)");
 544        kaonspectra->Draw();
 545
 546        c9->cd(4);
 547        chargedspectra->SetFillColor(42);
 548        chargedspectra->SetXTitle("(GeV)");
 549        chargedspectra->Draw();
 550
 551        break;
 552
 553      case 3:
 554
 555        if (nrawclusters) {
 556          TCanvas *c3=new TCanvas("c3","Clusters Statistics",50,50,600,700);
 557          c3->Divide(2,2);
 558
 559          c3->cd(1);
 560          c3->SetLogy(1);
 561          Clcharge->SetFillColor(42);
 562          Clcharge->SetXTitle("ADC units");
 563          Clcharge->Draw();
 564
 565          c3->cd(2);
 566          padnumber->SetFillColor(42);
 567          padnumber->SetXTitle("(counts)");
 568          padnumber->Draw();
 569
 570          c3->cd(3);
 571          clusev->SetFillColor(42);
 572          clusev->SetXTitle("(counts)");
 573          clusev->Draw();
 574
 575          c3->cd(4);
 576          padsmip->SetFillColor(42);
 577          padsmip->SetXTitle("(counts)");
 578          padsmip->Draw();
 579        }
 580
 581        if (nev<1)
 582          {
 583            TCanvas *c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
 584            mother->SetFillColor(42);
 585            mother->SetXTitle("counts");
 586            mother->Draw();
 587          }
 588
 589        TCanvas *c7 = new TCanvas("c7","Production Statistics",100,100,600,700);
 590        c7->Divide(2,2);
 591
 592        c7->cd(1);
 593        totalphotonsevent->SetFillColor(42);
 594        totalphotonsevent->SetXTitle("Photons (counts)");
 595        totalphotonsevent->Draw();
 596
 597        c7->cd(2);
 598        photev->SetFillColor(42);
 599        photev->SetXTitle("(counts)");
 600        photev->Draw();
 601
 602        c7->cd(3);
 603        feedev->SetFillColor(42);
 604        feedev->SetXTitle("(counts)");
 605        feedev->Draw();
 606
 607        c7->cd(4);
 608        padsev->SetFillColor(42);
 609        padsev->SetXTitle("(counts)");
 610        padsev->Draw();
 611
 612        break;
 613
 614      case 4:
 615
 616        TCanvas *c2 = new TCanvas("c2","Angles of incidence",50,50,600,700);
 617        c2->Divide(2,2);
 618
 619        c2->cd(1);
 620        hitsPhi->SetFillColor(42);
 621        hitsPhi->Draw();
 622        c2->cd(2);
 623        hitsTheta->SetFillColor(42);
 624        hitsTheta->Draw();
 625        c2->cd(3);
 626        Phi->SetFillColor(42);
 627        Phi->Draw();
 628        c2->cd(4);
 629        Theta->SetFillColor(42);
 630        Theta->Draw();
 631
 632
 633        TCanvas *c5 = new TCanvas("c5","Ring Reconstruction",100,100,900,700);
 634        c5->Divide(3,2);
 635
 636        c5->cd(1);
 637        ckovangle->SetFillColor(42);
 638        ckovangle->SetXTitle("angle (radians)");
 639        ckovangle->Draw();
 640
 641        c5->cd(2);
 642        radius->SetFillColor(42);
 643        radius->SetXTitle("radius (cm)");
 644        radius->Draw();
 645
 646        c5->cd(3);
 647        hc0->SetXTitle("pads");
 648        hc0->Draw("box");
 649
 650        c5->cd(5);
 651        Omega->SetFillColor(42);
 652        Omega->SetXTitle("angle (radians)");
 653        Omega->Draw();
 654
 655        c5->cd(4);
 656        PhotonCer->SetFillColor(42);
 657        PhotonCer->SetXTitle("angle (radians)");
 658        PhotonCer->Draw();
 659
 660        c5->cd(6);
 661        PadsUsed->SetXTitle("pads");
 662        PadsUsed->Draw("box");
 663
 664        break;
 665
 666      case 5:
 667
 668        if (ndigits)
 669          {
 670            TCanvas *c1 = new TCanvas("c1","Alice RICH digits",50,50,1200,700);
 671            c1->Divide(4,2);
 672            c1->cd(1);
 673            hc1->SetXTitle("ix (npads)");
 674            hc1->Draw("box");
 675            c1->cd(2);
 676            hc2->SetXTitle("ix (npads)");
 677            hc2->Draw("box");
 678            c1->cd(3);
 679            hc3->SetXTitle("ix (npads)");
 680            hc3->Draw("box");
 681            c1->cd(4);
 682            hc4->SetXTitle("ix (npads)");
 683            hc4->Draw("box");
 684            c1->cd(5);
 685            hc5->SetXTitle("ix (npads)");
 686            hc5->Draw("box");
 687            c1->cd(6);
 688            hc6->SetXTitle("ix (npads)");
 689            hc6->Draw("box");
 690            c1->cd(7);
 691            hc7->SetXTitle("ix (npads)");
 692            hc7->Draw("box");
 693            c1->cd(8);
 694            hc0->SetXTitle("ix (npads)");
 695            hc0->Draw("box");
 696          }
 697 //
 698        TCanvas *c4 = new TCanvas("c4","Hits per type",100,100,600,700);
 699        c4->Divide(2,2);
 700
 701        c4->cd(1);
 702        feedback->SetXTitle("x (cm)");
 703        feedback->SetYTitle("y (cm)");
 704        feedback->Draw();
 705
 706        c4->cd(2);
 707        //mip->SetFillColor(42);
 708        mip->SetXTitle("x (cm)");
 709        mip->SetYTitle("y (cm)");
 710        mip->Draw();
 711
 712        c4->cd(3);
 713        //cerenkov->SetFillColor(42);
 714        cerenkov->SetXTitle("x (cm)");
 715        cerenkov->SetYTitle("y (cm)");
 716        cerenkov->Draw();
 717
 718        c4->cd(4);
 719        //h->SetFillColor(42);
 720        h->SetXTitle("x (cm)");
 721        h->SetYTitle("y (cm)");
 722        h->Draw();
 723
 724        TCanvas *c10 = new TCanvas("c10","Hits distribution",150,150,600,350);
 725        c10->Divide(2,1);
 726
 727        c10->cd(1);
 728        hitsX->SetFillColor(42);
 729        hitsX->SetXTitle("(cm)");
 730        hitsX->Draw();
 731
 732        c10->cd(2);
 733        hitsY->SetFillColor(42);
 734        hitsY->SetXTitle("(cm)");
 735        hitsY->Draw();
 736
 737
 738        break;
 739
 740      }
 741
 742
 743    // calculate the number of pads which give a signal
 744
 745
 746    Int_t Np=0;
 747    for (Int_t i=0;i< NpadX;i++) {
 748        for (Int_t j=0;j< NpadY;j++) {
 749            if (Pad[i][j]>=6){
 750                Np+=1;
 751            }
 752        }
 753    }
 754    //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
 755    printf("End of macro\n");
 756 }
 757
 758