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