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