/////////////////////////////////////////////////////////////////////////
- Int_t NpadX = 162; // number of pads on X
- Int_t NpadY = 162; // number of pads on Y
-
- Int_t Pad[162][162];
- for (Int_t i=0;i<NpadX;i++) {
- for (Int_t j=0;j<NpadY;j++) {
- Pad[i][j]=0;
- }
- }
-
-
-
// Dynamically link some shared libs
if (gClassTable->GetID("AliRun") < 0) {
if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
}
-// Create some histograms
-
- Int_t xmin= -NpadX/2;
- Int_t xmax= NpadX/2;
- Int_t ymin= -NpadY/2;
- Int_t ymax= NpadY/2;
-
- TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
- TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
- TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
- TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
- TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
- TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
- TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
- TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
- TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
- TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
- TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
- TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
- TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
- TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
- TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
- TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
- TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
- TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
- TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
- TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
- TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
- TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
-
-
-
-
-// Start loop over events
-
- Int_t Nh=0;
- Int_t pads=0;
- Int_t Nh1=0;
- Float_t mom[3];
- Int_t nraw=0;
- Int_t phot=0;
- Int_t feed=0;
- Int_t padmip=0;
- Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
- Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
- Int_t chargedmuons=0, photons=0, primaryphotons=0, highprimaryphotons=0;
- TRandom random;
-
- for (int nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
-
-
- printf ("Event number : %d\n",nev);
- printf ("Number of particles: %d\n",nparticles);
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
-// Get pointers to RICH detector and Hits containers
+ // Get pointers to RICH detector and Hits containers
- AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
- Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
- gAlice->TreeR()->GetEvent(nent-1);
- TClonesArray *Rawclusters = RICH->RawClustAddress(2); // Raw clusters branch
- Int_t nrawclusters = Rawclusters->GetEntriesFast();
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
-
-
-
-// Start loop on tracks in the hits containers
- Int_t Nc=0;
- for (Int_t track=0; track<ntracks;track++) {
- printf ("Processing Track: %d\n",track);
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- if (RICH) {
- TClonesArray *PadHits = RICH->PadHits(); // Cluster branch
- TClonesArray *Hits = RICH->Hits(); // Hits branch
- TClonesArray *Cerenkovs = RICH->Cerenkovs(); // Cerenkovs branch
- }
- //see hits distribution
- Int_t nhits = Hits->GetEntriesFast();
- if (nhits) Nh+=nhits;
- for (Int_t hit=0;hit<nhits;hit++) {
- mHit = (AliRICHHit*) Hits->UncheckedAt(hit);
- Int_t nch = mHit->fChamber; // chamber number
- Float_t x = mHit->X(); // x-pos of hit
- Float_t y = mHit->Z(); // y-pos
- Float_t z = mHit->Y();
- Int_t index = mHit->Track();
- Int_t particle = mHit->fParticle;
- Float_t R;
-
- TParticle *current = (TParticle*)(*gAlice->Particles())[index];
-
- Float_t energy=current->Energy();
-
- R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
-
- //printf("Particle type: %d\n",current->GetPdgCode());
- if (TMath::Abs(particle) < 50000051)
- {
- //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
- if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
- {
- //gMC->Rndm(&random, 1);
- if (random->Rndm() < .1)
- production->Fill(current->Vz(),R,(float) 1);
- if (TMath::Abs(particle) == 50000050)
- {
- photons +=1;
- if (R<.005)
- {
- primaryphotons +=1;
- if (current->Energy()>0.001)
- highprimaryphotons +=1;
- }
- }
- if (TMath::Abs(particle) == 2112)
- {
- neutron +=1;
- if (current->Energy()>0.0001)
- highneutrons +=1;
- }
- }
- else
- {
- production->Fill(current->Vz(),R,(float) 1);
- //printf("Adding %d at %f\n",particle,R);
- }
- //mip->Fill(x,y,(float) 1);
- }
-
- if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
- {
- pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
- pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (R>2.5 && R<4.5)
- {
- pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- //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);
- }
- //printf("Pion mass: %e\n",current->GetCalcMass());
- pion +=1;
- if (TMath::Abs(particle)==211)
- {
- chargedpions +=1;
- if (R<.005)
- {
- primarypions +=1;
- if (current->Energy()>1)
- highprimarypions +=1;
- }
- }
- }
- if (TMath::Abs(particle)==2212)
- {
- protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
- protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (R>3 && R<4.3)
- protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass());
- proton +=1;
- }
- if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
- || TMath::Abs(particle)==311)
- {
- kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
- kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (R>2.5 && R<4.5)
- kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- //printf("Kaon mass: %e\n",current->GetCalcMass());
- kaon +=1;
- if (TMath::Abs(particle)==321)
- {
- chargedkaons +=1;
- if (R<.005)
- {
- primarykaons +=1;
- if (current->Energy()>1)
- highprimarykaons +=1;
- }
- }
- }
- if (TMath::Abs(particle)==11)
- {
- electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
- electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (R>2.5 && R<4.5)
- electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- //printf("Electron mass: %e\n",current->GetCalcMass());
- if (particle == 11)
- electron +=1;
- if (particle == -11)
- positron +=1;
- }
- if (TMath::Abs(particle)==13)
- {
- muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
- muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (R>2.5 && R<4.5)
- muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- //printf("Muon mass: %e\n",current->GetCalcMass());
- muon +=1;
- }
- if (TMath::Abs(particle)==2112)
- {
- neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
- neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (R>2.5 && R<4.5)
- {
- neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- //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);
- }
- //printf("Neutron mass: %e\n",current->GetCalcMass());
- neutron +=1;
- }
- if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
- {
- if (current->Energy()-current->GetCalcMass()>1)
- {
- chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
- chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- if (R>2.5 && R<4.5)
- chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
- }
- }
- //printf("Hits:%d\n",hit);
- //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
- // Fill the histograms
- Nh1+=nhits;
- //h->Fill(x,y,(float) 1);
- //}
- //}
- }
-
- }
-
- }
-
- //Create canvases, set the view range, show histograms
-
- TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
- production->SetFillColor(42);
- production->SetXTitle("z (m)");
- production->SetYTitle("R (m)");
- production->Draw();
-
- TCanvas *c16 = new TCanvas("c16","Particles Spectra II",100,100,600,350);
- c16->Divide(2,1);
-
- c16->cd(1);
- //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
- electronspectra1->SetFillColor(42);
- electronspectra1->SetXTitle("log(GeV)");
- electronspectra2->SetFillColor(46);
- electronspectra2->SetXTitle("log(GeV)");
- electronspectra3->SetFillColor(10);
- electronspectra3->SetXTitle("log(GeV)");
- //c13->SetLogx();
- electronspectra1->Draw();
- electronspectra2->Draw("same");
- electronspectra3->Draw("same");
-
- c16->cd(2);
- //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
- muonspectra1->SetFillColor(42);
- muonspectra1->SetXTitle("log(GeV)");
- muonspectra2->SetFillColor(46);
- muonspectra2->SetXTitle("log(GeV)");
- muonspectra3->SetFillColor(10);
- muonspectra3->SetXTitle("log(GeV)");
- //c14->SetLogx();
- muonspectra1->Draw();
- muonspectra2->Draw("same");
- muonspectra3->Draw("same");
-
- //c16->cd(3);
- //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
- //neutronspectra1->SetFillColor(42);
- //neutronspectra1->SetXTitle("log(GeV)");
- //neutronspectra2->SetFillColor(46);
- //neutronspectra2->SetXTitle("log(GeV)");
- //neutronspectra3->SetFillColor(10);
- //neutronspectra3->SetXTitle("log(GeV)");
- //c16->SetLogx();
- //neutronspectra1->Draw();
- //neutronspectra2->Draw("same");
- //neutronspectra3->Draw("same");
-
- TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
- //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
- c9->Divide(2,2);
-
- c9->cd(1);
- pionspectra1->SetFillColor(42);
- pionspectra1->SetXTitle("log(GeV)");
- pionspectra2->SetFillColor(46);
- pionspectra2->SetXTitle("log(GeV)");
- pionspectra3->SetFillColor(10);
- pionspectra3->SetXTitle("log(GeV)");
- //c9->SetLogx();
- pionspectra1->Draw();
- pionspectra2->Draw("same");
- pionspectra3->Draw("same");
-
- c9->cd(2);
- //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
- protonspectra1->SetFillColor(42);
- protonspectra1->SetXTitle("log(GeV)");
- protonspectra2->SetFillColor(46);
- protonspectra2->SetXTitle("log(GeV)");
- protonspectra3->SetFillColor(10);
- protonspectra3->SetXTitle("log(GeV)");
- //c10->SetLogx();
- protonspectra1->Draw();
- protonspectra2->Draw("same");
- protonspectra3->Draw("same");
-
- c9->cd(3);
- //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
- kaonspectra1->SetFillColor(42);
- kaonspectra1->SetXTitle("log(GeV)");
- kaonspectra2->SetFillColor(46);
- kaonspectra2->SetXTitle("log(GeV)");
- kaonspectra3->SetFillColor(10);
- kaonspectra3->SetXTitle("log(GeV)");
- //c11->SetLogx();
- kaonspectra1->Draw();
- kaonspectra2->Draw("same");
- kaonspectra3->Draw("same");
-
- c9->cd(4);
- //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
- chargedspectra1->SetFillColor(42);
- chargedspectra1->SetXTitle("log(GeV)");
- chargedspectra2->SetFillColor(46);
- chargedspectra2->SetXTitle("log(GeV)");
- chargedspectra3->SetFillColor(10);
- chargedspectra3->SetXTitle("log(GeV)");
- //c12->SetLogx();
- chargedspectra1->Draw();
- chargedspectra2->Draw("same");
- chargedspectra3->Draw("same");
-
- // calculate the number of pads which give a signal
-
-
- Int_t Np=0;
- for (Int_t i=0; i< NpadX;i++) {
- for (Int_t j=0;j< NpadY;j++) {
- if (Pad[i][j]>=6){
- Np+=1;
- }
- }
- }
- //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
-
- printf("*****************************************\n");
- printf("* Particle * Flux(m2) *\n");
- printf("*****************************************\n");
-
- printf("* Pions: * %3.1f *\n",pion/11.757312);
- printf("* Charged Pions: * %3.1f *\n",chargedpions/11.757312);
- printf("* Primary Pions: * %3.1f *\n",primarypions/11.757312);
- printf("* Primary Pions (p>1GeV/c): * %3.1f *\n",highprimarypions/11.757312);
- printf("* Kaons: * %3.1f *\n",kaon/11.757312);
- printf("* Charged Kaons: * %3.1f *\n",chargedkaons/11.757312);
- printf("* Primary Kaons: * %3.1f *\n",primarykaons/11.757312);
- printf("* Primary Kaons (p>1GeV/c): * %3.1f *\n",highprimarykaons/11.757312);
- printf("* Muons: * %3.1f *\n",muon/11.757312);
- printf("* Electrons: * %3.1f *\n",electron/11.757312);
- printf("* Positrons: * %3.1f *\n",positron/11.757312);
- printf("* Protons: * %3.1f *\n",proton/11.757312);
- printf("* All Charged: * %3.1f *\n",(chargedpions+chargedkaons+muon+electron+positron+proton)/11.757312);
- printf("*****************************************\n");
- //printf("* Photons: * %3.1f *\n",photons/11.757312);
- //printf("* Primary Photons: * %3.1f *\n",primaryphotons/11.757312);
- //printf("* Primary Photons (p>1MeV/c):* %3.1f *\n",highprimaryphotons/11.757312);
- //printf("*****************************************\n");
- //printf("* Neutrons: * %3.1f *\n",neutron);
- //printf("* Neutrons (p>100keV/c): * %3.1f *\n",highneutrons);
- //printf("*****************************************\n");
-
- printf("\nEnd of macro\n");
+ AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
+
+ RICH->DiagnosticsFE(evNumber1,evNumber2);
}
-
-
-
-
-