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