RICH/RICHpadtestC.C

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