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     else {
  33       //delete gAlice;
  34       gAlice = 0;
  35     }
  36
  37     gAlice=0;
  38
  39 // Connect the Root Galice file containing Geometry, Kine and Hits
  40
  41     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  42     if (!file) file = new TFile("galice.root","UPDATE");
  43
  44 // Get AliRun object from file or create it if not on file
  45
  46     if (!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     else {
  52       delete gAlice;
  53       gAlice = (AliRun*)file->Get("gAlice");
  54       if (gAlice) printf("AliRun object found on file\n");
  55       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  56     }
  57
  58 //  Create some histograms
  59
  60    Int_t xmin= -NpadX/2;
  61    Int_t xmax=  NpadX/2;
  62    Int_t ymin= -NpadY/2;
  63    Int_t ymax=  NpadY/2;
  64
  65    TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
  66    TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
  67    TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
  68    TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
  69    TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
  70    TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
  71    TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
  72    TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
  73    TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
  74    TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
  75    TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
  76    TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
  77    TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
  78    TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
  79    TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
  80    TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
  81    TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
  82    TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
  83    TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
  84    TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
  85    TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
  86    TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
  87
  88
  89
  90
  91 //   Start loop over events
  92
  93    Int_t Nh=0;
  94    Int_t pads=0;
  95    Int_t Nh1=0;
  96    Float_t mom[3];
  97    Int_t nraw=0;
  98    Int_t phot=0;
  99    Int_t feed=0;
 100    Int_t padmip=0;
 101    Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
 102    Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
 103    Int_t chargedmuons=0, photons=0, primaryphotons=0, highprimaryphotons=0;
 104
 105    TRandom random;
 106
 107    for (int nev=0; nev<= evNumber2; nev++) {
 108        Int_t nparticles = gAlice->GetEvent(nev);
 109
 110
 111        printf ("Event number       : %d\n",nev);
 112        printf ("Number of particles: %d\n",nparticles);
 113        if (nev < evNumber1) continue;
 114        if (nparticles <= 0) return;
 115
 116 // Get pointers to RICH detector and Hits containers
 117
 118        AliRICH *RICH  = (AliRICH*)gAlice->GetDetector("RICH");
 119        Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
 120        gAlice->TreeR()->GetEvent(nent-1);
 121        TClonesArray *Rawclusters = RICH->RawClustAddress(2);    //  Raw clusters branch
 122        Int_t nrawclusters = Rawclusters->GetEntriesFast();
 123        TTree *TH = gAlice->TreeH();
 124        Int_t ntracks = TH->GetEntries();
 125
 126
 127
 128 // Start loop on tracks in the hits containers
 129        Int_t Nc=0;
 130        for (Int_t track=0; track<ntracks;track++) {
 131            printf ("Processing Track: %d\n",track);
 132            gAlice->ResetHits();
 133            Int_t nbytes += TH->GetEvent(track);
 134            if (RICH)  {
 135                TClonesArray *PadHits = RICH->PadHits();      // Cluster branch
 136                TClonesArray *Hits = RICH->Hits();            // Hits branch
 137                TClonesArray *Cerenkovs = RICH->Cerenkovs();  // Cerenkovs branch
 138            }
 139            //see hits distribution
 140            Int_t nhits = Hits->GetEntriesFast();
 141            if (nhits) Nh+=nhits;
 142            for (Int_t hit=0;hit<nhits;hit++) {
 143               mHit = (AliRICHHit*) Hits->UncheckedAt(hit);
 144               Int_t nch  = mHit->fChamber;              // chamber number
 145               Float_t x  = mHit->X();                    // x-pos of hit
 146               Float_t y  = mHit->Z();                    // y-pos
 147               Float_t z  = mHit->Y();
 148               Int_t index = mHit->Track();
 149               Int_t particle = mHit->fParticle;
 150               Float_t R;
 151
 152               TParticle *current = (TParticle*)(*gAlice->Particles())[index];
 153
 154               Float_t energy=current->Energy();
 155
 156               R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
 157
 158               //printf("Particle type: %d\n",current->GetPdgCode());
 159               if (TMath::Abs(particle) < 50000051)
 160                 {
 161                   //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
 162                   if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
 163                     {
 164                       //gMC->Rndm(&random, 1);
 165                       if (random->Rndm() < .1)
 166                         production->Fill(current->Vz(),R,(float) 1);
 167                       if (TMath::Abs(particle) == 50000050)
 168                         {
 169                           photons +=1;
 170                           if (R<.005)
 171                             {
 172                               primaryphotons +=1;
 173                               if (current->Energy()>0.001)
 174                                 highprimaryphotons +=1;
 175                             }
 176                         }
 177                       if (TMath::Abs(particle) == 2112)
 178                         {
 179                           neutron +=1;
 180                           if (current->Energy()>0.0001)
 181                             highneutrons +=1;
 182                         }
 183                     }
 184                   else
 185                     {
 186                       production->Fill(current->Vz(),R,(float) 1);
 187                       //printf("Adding %d at %f\n",particle,R);
 188                     }
 189                   //mip->Fill(x,y,(float) 1);
 190                 }
 191
 192               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
 193                 {
 194                   pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 195                   if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
 196                     pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 197                   if (R>2.5 && R<4.5)
 198                     {
 199                     pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 200                     //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);
 201                     }
 202                   //printf("Pion mass: %e\n",current->GetCalcMass());
 203                   pion +=1;
 204                   if (TMath::Abs(particle)==211)
 205                     {
 206                       chargedpions +=1;
 207                       if (R<.005)
 208                         {
 209                           primarypions +=1;
 210                           if (current->Energy()>1)
 211                             highprimarypions +=1;
 212                         }
 213                     }
 214                 }
 215               if (TMath::Abs(particle)==2212)
 216                 {
 217                   protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 218                   if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
 219                     protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 220                   if (R>3 && R<4.3)
 221                     protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 222                   //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass());
 223                   proton +=1;
 224                 }
 225               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
 226                   || TMath::Abs(particle)==311)
 227                 {
 228                   kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 229                   if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
 230                     kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 231                   if (R>2.5 && R<4.5)
 232                     kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 233                   //printf("Kaon mass: %e\n",current->GetCalcMass());
 234                   kaon +=1;
 235                   if (TMath::Abs(particle)==321)
 236                     {
 237                       chargedkaons +=1;
 238                       if (R<.005)
 239                         {
 240                           primarykaons +=1;
 241                           if (current->Energy()>1)
 242                             highprimarykaons +=1;
 243                         }
 244                     }
 245                 }
 246               if (TMath::Abs(particle)==11)
 247                 {
 248                   electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 249                   if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
 250                     electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 251                   if (R>2.5 && R<4.5)
 252                     electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 253                   //printf("Electron mass: %e\n",current->GetCalcMass());
 254                   if (particle == 11)
 255                     electron +=1;
 256                   if (particle == -11)
 257                     positron +=1;
 258                 }
 259               if (TMath::Abs(particle)==13)
 260                 {
 261                   muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 262                   if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
 263                     muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 264                   if (R>2.5 && R<4.5)
 265                     muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 266                   //printf("Muon mass: %e\n",current->GetCalcMass());
 267                   muon +=1;
 268                 }
 269               if (TMath::Abs(particle)==2112)
 270                 {
 271                   neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 272                   if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
 273                     neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 274                   if (R>2.5 && R<4.5)
 275                     {
 276                       neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 277                       //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);
 278                     }
 279                   //printf("Neutron mass: %e\n",current->GetCalcMass());
 280                   neutron +=1;
 281                 }
 282               if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
 283                 {
 284                   if (current->Energy()-current->GetCalcMass()>1)
 285                     {
 286                       chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 287                       if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
 288                         chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 289                       if (R>2.5 && R<4.5)
 290                         chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
 291                     }
 292                 }
 293               //printf("Hits:%d\n",hit);
 294               //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
 295               // Fill the histograms
 296               Nh1+=nhits;
 297               //h->Fill(x,y,(float) 1);
 298                   //}
 299               //}
 300            }
 301
 302        }
 303
 304    }
 305
 306    //Create canvases, set the view range, show histograms
 307
 308    TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
 309    production->SetFillColor(42);
 310    production->SetXTitle("z (m)");
 311    production->SetYTitle("R (m)");
 312    production->Draw();
 313
 314    TCanvas *c16 = new TCanvas("c16","Particles Spectra II",100,100,600,350);
 315    c16->Divide(2,1);
 316
 317    c16->cd(1);
 318    //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
 319    electronspectra1->SetFillColor(42);
 320    electronspectra1->SetXTitle("log(GeV)");
 321    electronspectra2->SetFillColor(46);
 322    electronspectra2->SetXTitle("log(GeV)");
 323    electronspectra3->SetFillColor(10);
 324    electronspectra3->SetXTitle("log(GeV)");
 325    //c13->SetLogx();
 326    electronspectra1->Draw();
 327    electronspectra2->Draw("same");
 328    electronspectra3->Draw("same");
 329
 330    c16->cd(2);
 331    //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
 332    muonspectra1->SetFillColor(42);
 333    muonspectra1->SetXTitle("log(GeV)");
 334    muonspectra2->SetFillColor(46);
 335    muonspectra2->SetXTitle("log(GeV)");
 336    muonspectra3->SetFillColor(10);
 337    muonspectra3->SetXTitle("log(GeV)");
 338    //c14->SetLogx();
 339    muonspectra1->Draw();
 340    muonspectra2->Draw("same");
 341    muonspectra3->Draw("same");
 342
 343    //c16->cd(3);
 344    //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
 345    //neutronspectra1->SetFillColor(42);
 346    //neutronspectra1->SetXTitle("log(GeV)");
 347    //neutronspectra2->SetFillColor(46);
 348    //neutronspectra2->SetXTitle("log(GeV)");
 349    //neutronspectra3->SetFillColor(10);
 350    //neutronspectra3->SetXTitle("log(GeV)");
 351    //c16->SetLogx();
 352    //neutronspectra1->Draw();
 353    //neutronspectra2->Draw("same");
 354    //neutronspectra3->Draw("same");
 355
 356    TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
 357    //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
 358    c9->Divide(2,2);
 359
 360    c9->cd(1);
 361    pionspectra1->SetFillColor(42);
 362    pionspectra1->SetXTitle("log(GeV)");
 363    pionspectra2->SetFillColor(46);
 364    pionspectra2->SetXTitle("log(GeV)");
 365    pionspectra3->SetFillColor(10);
 366    pionspectra3->SetXTitle("log(GeV)");
 367    //c9->SetLogx();
 368    pionspectra1->Draw();
 369    pionspectra2->Draw("same");
 370    pionspectra3->Draw("same");
 371
 372    c9->cd(2);
 373    //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
 374    protonspectra1->SetFillColor(42);
 375    protonspectra1->SetXTitle("log(GeV)");
 376    protonspectra2->SetFillColor(46);
 377    protonspectra2->SetXTitle("log(GeV)");
 378    protonspectra3->SetFillColor(10);
 379    protonspectra3->SetXTitle("log(GeV)");
 380    //c10->SetLogx();
 381    protonspectra1->Draw();
 382    protonspectra2->Draw("same");
 383    protonspectra3->Draw("same");
 384
 385    c9->cd(3);
 386    //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
 387    kaonspectra1->SetFillColor(42);
 388    kaonspectra1->SetXTitle("log(GeV)");
 389    kaonspectra2->SetFillColor(46);
 390    kaonspectra2->SetXTitle("log(GeV)");
 391    kaonspectra3->SetFillColor(10);
 392    kaonspectra3->SetXTitle("log(GeV)");
 393    //c11->SetLogx();
 394    kaonspectra1->Draw();
 395    kaonspectra2->Draw("same");
 396    kaonspectra3->Draw("same");
 397
 398    c9->cd(4);
 399    //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
 400    chargedspectra1->SetFillColor(42);
 401    chargedspectra1->SetXTitle("log(GeV)");
 402    chargedspectra2->SetFillColor(46);
 403    chargedspectra2->SetXTitle("log(GeV)");
 404    chargedspectra3->SetFillColor(10);
 405    chargedspectra3->SetXTitle("log(GeV)");
 406    //c12->SetLogx();
 407    chargedspectra1->Draw();
 408    chargedspectra2->Draw("same");
 409    chargedspectra3->Draw("same");
 410
 411    // calculate the number of pads which give a signal
 412
 413
 414    Int_t Np=0;
 415    for (Int_t i=0; i< NpadX;i++) {
 416        for (Int_t j=0;j< NpadY;j++) {
 417            if (Pad[i][j]>=6){
 418                Np+=1;
 419            }
 420        }
 421    }
 422    //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
 423
 424    printf("*****************************************\n");
 425    printf("* Particle                   * Flux(m2) *\n");
 426    printf("*****************************************\n");
 427
 428    printf("* Pions:                     *   %3.1f   *\n",pion/11.757312);
 429    printf("* Charged Pions:             *   %3.1f   *\n",chargedpions/11.757312);
 430    printf("* Primary Pions:             *   %3.1f   *\n",primarypions/11.757312);
 431    printf("* Primary Pions (p>1GeV/c):  *   %3.1f   *\n",highprimarypions/11.757312);
 432    printf("* Kaons:                     *   %3.1f   *\n",kaon/11.757312);
 433    printf("* Charged Kaons:             *   %3.1f   *\n",chargedkaons/11.757312);
 434    printf("* Primary Kaons:             *   %3.1f   *\n",primarykaons/11.757312);
 435    printf("* Primary Kaons (p>1GeV/c):  *   %3.1f   *\n",highprimarykaons/11.757312);
 436    printf("* Muons:                     *   %3.1f   *\n",muon/11.757312);
 437    printf("* Electrons:                 *   %3.1f   *\n",electron/11.757312);
 438    printf("* Positrons:                 *   %3.1f   *\n",positron/11.757312);
 439    printf("* Protons:                   *   %3.1f   *\n",proton/11.757312);
 440    printf("* All Charged:               *   %3.1f   *\n",(chargedpions+chargedkaons+muon+electron+positron+proton)/11.757312);
 441    printf("*****************************************\n");
 442    //printf("* Photons:                   *   %3.1f   *\n",photons/11.757312);
 443    //printf("* Primary Photons:           *   %3.1f   *\n",primaryphotons/11.757312);
 444    //printf("* Primary Photons (p>1MeV/c):*   %3.1f   *\n",highprimaryphotons/11.757312);
 445    //printf("*****************************************\n");
 446    //printf("* Neutrons:                  *   %3.1f   *\n",neutron);
 447    //printf("* Neutrons (p>100keV/c):     *   %3.1f   *\n",highneutrons);
 448    //printf("*****************************************\n");
 449
 450    printf("\nEnd of macro\n");
 451 }
 452
 453
 454
 455
 456