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