[u/mrichter/AliRoot.git] / MUON / MUONhits.C

void MUONhits (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ic=1) 
{
/////////////////////////////////////////////////////////////////////////
//   This macro is a small example of a ROOT macro
//   illustrating how to read the output of GALICE
//   and do some analysis.
//   
/////////////////////////////////////////////////////////////////////////

// Dynamically link some shared libs                    
    Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
			56.1, 56.1, 59.6, 59.6};
    Float_t rmax[14] = {81.6, 81.6, 109.3, 109.3, 154.4, 154.4, 197.8, 
			197.8, 229.5, 229.5, 254., 254, 270., 270.};
    

    if (gClassTable->GetID("AliRun") < 0) {
	gROOT->LoadMacro("loadlibs.C");
	loadlibs();
    }
// Connect the Root Galice file containing Geometry, Kine and Hits

    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");

    
    if (!file) {
	printf("\n Creating galice.root \n");
	file = new TFile("galice.root");
    } else {
	printf("\n galice.root found in file list");
    }
    file->ls();
//
    TDatabasePDG*  DataBase = new TDatabasePDG();
    
// Get AliRun object from file or create it if not on file
    if (!gAlice) {
	gAlice = (AliRun*)(file->Get("gAlice"));
	if (gAlice) printf("AliRun object found on file\n");
	if (!gAlice) {
	    printf("\n create new gAlice object");
	    gAlice = new AliRun("gAlice","Alice test program");
	}
    }
    
//  Create some histograms

    TH1F *Hits[14],   *HitDensity[14];
    TH1F *Hits_g[14], *HitDensity_g[14];
    TH1F *Hits_n[14], *HitDensity_n[14];
    TH1F *Mom[14], *Mom_g[14], *Mom_n[14];
    
    for (Int_t i=0; i<14; i++) {
	Hits[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
	Hits[i]->SetFillColor(0);
	Hits[i]->SetXTitle("R (cm)");
	
	HitDensity[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
	HitDensity[i]->SetFillColor(0);
	HitDensity[i]->SetXTitle("R (cm)");

	Hits_g[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
	Hits_g[i]->SetFillColor(0);
	Hits_g[i]->SetXTitle("R (cm)");

	HitDensity_g[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
	HitDensity_g[i]->SetFillColor(0);
	HitDensity_g[i]->SetXTitle("R (cm)");

	Mom[i]  =  new TH1F("mom","Energy",70,-6,1);
	Mom[i]->SetFillColor(0);
	Mom[i]->SetXTitle("E (GeV)");

	Mom_g[i]  =  new TH1F("mom","Energy",70,-6,1);
	Mom_g[i]->SetFillColor(0);
	Mom_g[i]->SetXTitle("E (GeV)");

	Mom_n[i]  =  new TH1F("mom","Energy",70,-6,1);
	Mom_n[i]->SetFillColor(0);
	Mom_n[i]->SetXTitle("E (GeV)");

	Hits_n[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
	Hits_n[i]->SetFillColor(0);
	Hits_n[i]->SetXTitle("R (cm)");

	HitDensity_n[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
	HitDensity_n[i]->SetFillColor(0);
	HitDensity_n[i]->SetXTitle("R (cm)");
    }

    TH1F *theta   =  new TH1F("theta","Theta distribution",180,0,180);
 
    TH1F *emult   =  new TH1F("emult","Event Multiplicity",100,0,1000);   
    AliMUONChamber*  iChamber;
    AliMUONSegmentation* seg;
    
//
//   Loop over events 
//
    Float_t dpx[2], dpy[2];
    Int_t mdig =0;
    Int_t nhit=0;
    
    for (Int_t nev=0; nev<= evNumber2; nev++) {
	Int_t nparticles = gAlice->GetEvent(nev);
	if (nev < evNumber1) continue;
	if (nparticles <= 0) return;
	
	AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
	printf("\n track %d %d \n ", nev, MUON);
	
	TTree *TH = gAlice->TreeH();
	Int_t ntracks = TH->GetEntries();
//
//   Loop over events
//
	Int_t EvMult=0;
       
	for (Int_t track=0; track<ntracks;track++) {
	    gAlice->ResetHits();
	    Int_t nbytes += TH->GetEvent(track);
	    if (MUON)  {
//
//  Loop over hits
//
		Int_t NperCh=0;
	       
		for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
		    mHit;
		    mHit=(AliMUONHit*)MUON->NextHit()) 
		{
		    Int_t   nch   =    mHit->fChamber;  // chamber number
		    Float_t x     =    mHit->fX;        // x-pos of hit
		    Float_t y     =    mHit->fY;        // y-pos
		    Float_t Eloss =    mHit->fEloss;
		    Float_t Theta =    mHit->fTheta;
		    Float_t Particle = mHit->fParticle;
		    Float_t P    =     
			TMath::Sqrt(mHit->fCxHit*mHit->fCxHit+
				    mHit->fCyHit*mHit->fCyHit+
				    mHit->fCzHit*mHit->fCzHit);
		    TParticlePDG* Part = DataBase->GetParticle(Particle);
		    Double_t mass = Part->Mass();
		    
		    if (nch >13) continue;
		    if (nch ==1) EvMult++;
		    if (mHit->fAge > 5.e-6) continue;
		   
		    Float_t r=TMath::Sqrt(x*x+y*y);
		    if (nch ==0) continue;
		    
		    if (r < rmin[nch-1]) continue;
		    if (r > rmax[nch-1]) continue;
		   
		    Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);

		    Float_t Ekin=TMath::Sqrt(P*P+mass*mass)-mass;
		    if (Particle == 2112) {
			Hits_n[nch-1]->Fill(r,(float) 1);
			HitDensity_n[nch-1]->Fill(r,wgt);
			Mom_n[nch-1]->Fill(TMath::Log10(Ekin),1);
		    } else if (Particle == 22) {
			Hits_g[nch-1]->Fill(r,(float) 1);
			HitDensity_g[nch-1]->Fill(r,wgt);
			Mom_g[nch-1]->Fill(TMath::Log10(Ekin),1);
		    } else {
			Hits[nch-1]->Fill(r,(float) 1);
			HitDensity[nch-1]->Fill(r,wgt);
			Mom[nch-1]->Fill(TMath::Log10(Ekin),1);
		    }
		} // hit loop
	    } // if MUON
	} // track loop
    }
    
//Create a canvas, set the view range, show histograms
    Int_t k;
    TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
    c1->Divide(2,4);
    for (k=0; k<7; k++) {
	c1->cd(k+1);
	HitDensity[2*k]->Draw();
    }

    TCanvas *c2 = new TCanvas("c2","Hit Densities (gamma)",400,10,600,700);
    c2->Divide(2,4);
    for (k=0; k<7; k++) {
	c2->cd(k+1);
	HitDensity_g[2*k]->Draw();
    }
    
    TCanvas *c3 = new TCanvas("c3","Hit Densities (neutron)",400,10,600,700);
    c3->Divide(2,4);
    for (k=0; k<7; k++) {
	c3->cd(k+1);
	HitDensity_n[2*k]->Draw();
    }

    TCanvas *c4 = new TCanvas("c4","Energy (charged)",400,10,600,700);
    c4->Divide(2,4);
    for (k=0; k<7; k++) {
	c4->cd(k+1);
	Mom[2*k]->Draw();
    }

    TCanvas *c5 = new TCanvas("c5","Energy (gamma)",400,10,600,700);
    c5->Divide(2,4);
    for (k=0; k<7; k++) {
	c5->cd(k+1);
	Mom_g[2*k]->Draw();
    }

    TCanvas *c6 = new TCanvas("c6","Energy (gamma)",400,10,600,700);
    c6->Divide(2,4);
    for (k=0; k<7; k++) {
	c6->cd(k+1);
	Mom_n[2*k]->Draw();
    }
}
Commit	Line	Data
a9e2aefa	1	void MUONhits (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ic=1)
	2	{
	3	/////////////////////////////////////////////////////////////////////////
	4	// This macro is a small example of a ROOT macro
	5	// illustrating how to read the output of GALICE
	6	// and do some analysis.
	7	//
	8	/////////////////////////////////////////////////////////////////////////
	9
	10	// Dynamically link some shared libs
	11	Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
	12	56.1, 56.1, 59.6, 59.6};
	13	Float_t rmax[14] = {81.6, 81.6, 109.3, 109.3, 154.4, 154.4, 197.8,
	14	197.8, 229.5, 229.5, 254., 254, 270., 270.};
	15
	16
	17	if (gClassTable->GetID("AliRun") < 0) {
	18	gROOT->LoadMacro("loadlibs.C");
	19	loadlibs();
	20	}
	21	// Connect the Root Galice file containing Geometry, Kine and Hits
	22
	23	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject("galice.root");
	24
	25
	26	if (!file) {
	27	printf("\n Creating galice.root \n");
	28	file = new TFile("galice.root");
	29	} else {
	30	printf("\n galice.root found in file list");
	31	}
	32	file->ls();
	33	//
	34	TDatabasePDG* DataBase = new TDatabasePDG();
	35
	36	// Get AliRun object from file or create it if not on file
	37	if (!gAlice) {
	38	gAlice = (AliRun*)(file->Get("gAlice"));
	39	if (gAlice) printf("AliRun object found on file\n");
	40	if (!gAlice) {
	41	printf("\n create new gAlice object");
	42	gAlice = new AliRun("gAlice","Alice test program");
	43	}
	44	}
	45
	46	// Create some histograms
	47
	48	TH1F Hits[14], HitDensity[14];
	49	TH1F Hits_g[14], HitDensity_g[14];
	50	TH1F Hits_n[14], HitDensity_n[14];
	51	TH1F Mom[14], Mom_g[14], *Mom_n[14];
	52
	53	for (Int_t i=0; i<14; i++) {
	54	Hits[i] = new TH1F("hits1","Hit Distribution",30,0,300);
	55	Hits[i]->SetFillColor(0);
	56	Hits[i]->SetXTitle("R (cm)");
	57
	58	HitDensity[i] = new TH1F("dhits1","Hit Density",30,0,300);
	59	HitDensity[i]->SetFillColor(0);
	60	HitDensity[i]->SetXTitle("R (cm)");
	61
	62	Hits_g[i] = new TH1F("hits1","Hit Distribution",30,0,300);
	63	Hits_g[i]->SetFillColor(0);
	64	Hits_g[i]->SetXTitle("R (cm)");
65
66	HitDensity_g[i] = new TH1F("dhits1","Hit Density",30,0,300);
67	HitDensity_g[i]->SetFillColor(0);
68	HitDensity_g[i]->SetXTitle("R (cm)");
69
70	Mom[i] = new TH1F("mom","Energy",70,-6,1);
71	Mom[i]->SetFillColor(0);
72	Mom[i]->SetXTitle("E (GeV)");
73
74	Mom_g[i] = new TH1F("mom","Energy",70,-6,1);
75	Mom_g[i]->SetFillColor(0);
76	Mom_g[i]->SetXTitle("E (GeV)");
77
78	Mom_n[i] = new TH1F("mom","Energy",70,-6,1);
79	Mom_n[i]->SetFillColor(0);
80	Mom_n[i]->SetXTitle("E (GeV)");
81
82	Hits_n[i] = new TH1F("hits1","Hit Distribution",30,0,300);
83	Hits_n[i]->SetFillColor(0);
84	Hits_n[i]->SetXTitle("R (cm)");
85
86	HitDensity_n[i] = new TH1F("dhits1","Hit Density",30,0,300);
87	HitDensity_n[i]->SetFillColor(0);
88	HitDensity_n[i]->SetXTitle("R (cm)");
89	}
90
91	TH1F *theta = new TH1F("theta","Theta distribution",180,0,180);
92
93	TH1F *emult = new TH1F("emult","Event Multiplicity",100,0,1000);
94	AliMUONChamber* iChamber;
95	AliMUONSegmentation* seg;
96
97	//
98	// Loop over events
99	//
100	Float_t dpx[2], dpy[2];
101	Int_t mdig =0;
102	Int_t nhit=0;
103
104	for (Int_t nev=0; nev<= evNumber2; nev++) {
105	Int_t nparticles = gAlice->GetEvent(nev);
106	if (nev < evNumber1) continue;
107	if (nparticles <= 0) return;
108
109	AliMUON MUON = (AliMUON) gAlice->GetModule("MUON");
110	printf("\n track %d %d \n ", nev, MUON);
111
112	TTree *TH = gAlice->TreeH();
113	Int_t ntracks = TH->GetEntries();
114	//
115	// Loop over events
116	//
117	Int_t EvMult=0;
118
119	for (Int_t track=0; track<ntracks;track++) {
120	gAlice->ResetHits();
121	Int_t nbytes += TH->GetEvent(track);
122	if (MUON) {
123	//
124	// Loop over hits
125	//
126	Int_t NperCh=0;
127
128	for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
129	mHit;
130	mHit=(AliMUONHit*)MUON->NextHit())
131	{
132	Int_t nch = mHit->fChamber; // chamber number
133	Float_t x = mHit->fX; // x-pos of hit
134	Float_t y = mHit->fY; // y-pos
135	Float_t Eloss = mHit->fEloss;
136	Float_t Theta = mHit->fTheta;
137	Float_t Particle = mHit->fParticle;
138	Float_t P =
139	TMath::Sqrt(mHit->fCxHit*mHit->fCxHit+
140	mHit->fCyHit*mHit->fCyHit+
141	mHit->fCzHit*mHit->fCzHit);
142	TParticlePDG* Part = DataBase->GetParticle(Particle);
143	Double_t mass = Part->Mass();
144
145	if (nch >13) continue;
146	if (nch ==1) EvMult++;
147	if (mHit->fAge > 5.e-6) continue;
148
149	Float_t r=TMath::Sqrt(xx+yy);
150	if (nch ==0) continue;
151
152	if (r < rmin[nch-1]) continue;
153	if (r > rmax[nch-1]) continue;
154
155	Float_t wgt=1/(210TMath::Pi()*r)/(evNumber2+1);
156
157	Float_t Ekin=TMath::Sqrt(PP+massmass)-mass;
158	if (Particle == 2112) {
159	Hits_n[nch-1]->Fill(r,(float) 1);
160	HitDensity_n[nch-1]->Fill(r,wgt);
161	Mom_n[nch-1]->Fill(TMath::Log10(Ekin),1);
162	} else if (Particle == 22) {
163	Hits_g[nch-1]->Fill(r,(float) 1);
164	HitDensity_g[nch-1]->Fill(r,wgt);
165	Mom_g[nch-1]->Fill(TMath::Log10(Ekin),1);
166	} else {
167	Hits[nch-1]->Fill(r,(float) 1);
168	HitDensity[nch-1]->Fill(r,wgt);
169	Mom[nch-1]->Fill(TMath::Log10(Ekin),1);
170	}
171	} // hit loop
172	} // if MUON
173	} // track loop
174	}
175
176	//Create a canvas, set the view range, show histograms
177	Int_t k;
178	TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
179	c1->Divide(2,4);
180	for (k=0; k<7; k++) {
181	c1->cd(k+1);
182	HitDensity[2*k]->Draw();
183	}
184
185	TCanvas *c2 = new TCanvas("c2","Hit Densities (gamma)",400,10,600,700);
186	c2->Divide(2,4);
187	for (k=0; k<7; k++) {
188	c2->cd(k+1);
189	HitDensity_g[2*k]->Draw();
190	}
191
192	TCanvas *c3 = new TCanvas("c3","Hit Densities (neutron)",400,10,600,700);
193	c3->Divide(2,4);
194	for (k=0; k<7; k++) {
195	c3->cd(k+1);
196	HitDensity_n[2*k]->Draw();
197	}
198
199	TCanvas *c4 = new TCanvas("c4","Energy (charged)",400,10,600,700);
200	c4->Divide(2,4);
201	for (k=0; k<7; k++) {
202	c4->cd(k+1);
203	Mom[2*k]->Draw();
204	}
205
206	TCanvas *c5 = new TCanvas("c5","Energy (gamma)",400,10,600,700);
207	c5->Divide(2,4);
208	for (k=0; k<7; k++) {
209	c5->cd(k+1);
210	Mom_g[2*k]->Draw();
211	}
212
213	TCanvas *c6 = new TCanvas("c6","Energy (gamma)",400,10,600,700);
214	c6->Divide(2,4);
215	for (k=0; k<7; k++) {
216	c6->cd(k+1);
217	Mom_n[2*k]->Draw();
218	}
219	}
220
221
222
223
224
225
226