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