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

void MUONocc (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                    

   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();
    
// 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[10], *HitDensity[10], *Occ0[10],*Occ1[10], *Mult[10];
    for (Int_t i=0; i<10; i++) {
	Hits[i]     =  new TH1F("hits1","Hit Distribution",26,0,260);
	Hits[i]->SetFillColor(0);
	Hits[i]->SetXTitle("R (cm)");

	HitDensity[i]  =  new TH1F("dhits1","Hit Density Distribution",26,0,260);
	HitDensity[i]->SetFillColor(0);
	HitDensity[i]->SetXTitle("R (cm)");

	Occ0[i]  =  new TH1F("occ0","Occupancy Density",26,0,260);
	Occ0[i] -> SetFillColor(0);
	Occ0[i] -> SetXTitle("R (cm)");
	Occ1[i]  =  new TH1F("occ1","Occupancy",26,0,260);
	Occ1[i] -> SetFillColor(0);
	Occ1[i] -> SetXTitle("R (cm)");
	Occ1[i] -> SetYTitle("Occupancy");
	Mult[i]  =  new TH1F("mult","Mult distribution",26,0,260);
	Mult[i] -> SetFillColor(0);
	Mult[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->Chamber();  // chamber number
		   Float_t x     = mHit->X();        // x-pos of hit
		   Float_t y     = mHit->Y();        // y-pos
		   Float_t Eloss = mHit->Eloss();
		   Float_t Theta = mHit->Theta();
		   theta->Fill(Theta,(float) 1);
		   if (nch >10) continue;
		   if (nch ==1) EvMult++;
//
//
		   iChamber = & MUON->Chamber(nch-1);
		   response= iChamber.GetResponseModel();
		   seg     = iChamber.GetSegmentationModel(1);
		   NperCh++;
		   Int_t i,j;
		   seg->GetPadIxy(x,y,i,j);
		   Int_t isec = seg->Sector(i,j);
		   if (isec ==1 && nch==ic) nhit++;
		   
		   Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
		   Float_t r=TMath::Sqrt(x*x+y*y);
		   Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);
		   wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
		   Hits[nch-1]->Fill(r,(float) 1);
		   HitDensity[nch-1]->Fill(r,wgt);
		   Occ0[nch-1]->Fill(r,wgt*a);
	       } // hit loop
	   } // if MUON
       } // track loop
       emult->Fill(Float_t(EvMult), (Float_t) 1);
       
       Int_t iseg=1;

       for (Int_t ich=0; ich<10; ich++) {
	   iChamber = & MUON->Chamber(ich);
	   seg=iChamber.GetSegmentationModel(iseg);
	   gAlice->ResetDigits();
	   gAlice->TreeD()->GetEvent(iseg);
	   TClonesArray *MUONDigits  = MUON->DigitsAddress(ich);
	   Int_t Ndigits=MUONDigits->GetEntriesFast();
	   AliMUONDigit* dig;
	   printf("\n Reading %d digits\n", Ndigits);
	   if (MUONDigits)  {
	       for (Int_t ndig=0; ndig<Ndigits; ndig++) 
	       {
		   dig = (AliMUONDigit*)MUONDigits->UncheckedAt(ndig); 
		   Int_t i=dig->PadX();
		   Int_t j=dig->PadY();
		   Float_t x,y;
		   seg->GetPadCxy(i,j,x,y);
		   Int_t isec = seg->Sector(i,j);
		   Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
		   Float_t r=TMath::Sqrt(x*x+y*y);
//	       if (r<25) 
//	       printf("\n Sector,a %d %f", isec,a);
		   if (isec==1) mdig++;
		   
		   Float_t wgt;
		   wgt=1/(2.*10.*TMath::Pi()*r)/(evNumber2+1)*a;
// Take into account inefficiency due to frames
		   wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
		   Occ1[ich]->Fill(r,wgt);
	       } // digit loop
	       Mult[ich]->Divide(Occ1[ich],Occ0[ich]);
	   } // chamber loop
       } // if MUONDigits
   } // event loop 
//
   printf("\n hits, digits %d %d\n ", nhit, mdig);
   

//Create a canvas, set the view range, show histograms
   TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
   pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
   pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
   pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
   pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
   pad11->SetFillColor(0);
   pad12->SetFillColor(0);
   pad13->SetFillColor(0);
   pad14->SetFillColor(0);
   pad11->Draw();
   pad12->Draw();
   pad13->Draw();
   pad14->Draw();
   
   pad11->cd();
   HitDensity[0]->Draw();
   pad12->cd();
   HitDensity[1]->Draw();
   pad13->cd();
   HitDensity[2]->Draw();
   pad14->cd();
   HitDensity[3]->Draw();

   TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
   pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
   pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
   pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
   pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
   pad21->SetFillColor(0);
   pad22->SetFillColor(0);
   pad23->SetFillColor(0);
   pad24->SetFillColor(0);
   pad21->Draw();
   pad22->Draw();
   pad23->Draw();
   pad24->Draw();

   pad21->cd();
   HitDensity[4]->Draw();
   pad22->cd();
   HitDensity[5]->Draw();
   pad23->cd();
   HitDensity[6]->Draw();
   pad24->cd();
   HitDensity[7]->Draw();

   TCanvas *c3 = new TCanvas("c3","Hit Densities",400,10,600,700);
   pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
   pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
   pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
   pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
   pad31->SetFillColor(0);
   pad32->SetFillColor(0);
   pad33->SetFillColor(0);
   pad34->SetFillColor(0);
   pad31->Draw();
   pad32->Draw();
   pad33->Draw();
   pad34->Draw();

   pad31->cd();
   HitDensity[8]->Draw();
   pad32->cd();
   HitDensity[9]->Draw();

   TCanvas *c4 = new TCanvas("c4","Occupancies",400,10,600,700);
   pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
   pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
   pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
   pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
   pad41->SetFillColor(0);
   pad42->SetFillColor(0);
   pad43->SetFillColor(0);
   pad44->SetFillColor(0);
   pad41->Draw();
   pad42->Draw();
   pad43->Draw();
   pad44->Draw();


   pad41->cd();
   Occ1[0]->Draw();
   pad42->cd();
   Occ1[2]->Draw();
   pad43->cd();
   Occ1[4]->Draw();
   pad44->cd();
   Occ1[6]->Scale(1.25);
   Occ1[6]->Draw();

   TCanvas *c5 = new TCanvas("c5","Occupancies",400,10,600,700);
   pad51 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
   pad52 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
   pad53 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
   pad54 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
   pad51->SetFillColor(0);
   pad52->SetFillColor(0);
   pad53->SetFillColor(0);
   pad54->SetFillColor(0);
   pad51->Draw();
   pad52->Draw();
   pad53->Draw();
   pad54->Draw();


   pad51->cd();
   Occ1[8]->Scale(1.25);
   
   Occ1[8]->Draw();
}
Commit	Line	Data
a9e2aefa	1	void MUONocc (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
	12	if (gClassTable->GetID("AliRun") < 0) {
	13	gROOT->LoadMacro("loadlibs.C");
	14	loadlibs();
	15	}
	16	// Connect the Root Galice file containing Geometry, Kine and Hits
	17
	18	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject("galice.root");
	19
	20
	21	if (!file) {
	22	printf("\n Creating galice.root \n");
	23	file = new TFile("galice.root");
	24	} else {
	25	printf("\n galice.root found in file list");
	26	}
	27	file->ls();
	28
	29	// Get AliRun object from file or create it if not on file
	30	if (!gAlice) {
	31	gAlice = (AliRun*)(file->Get("gAlice"));
	32	if (gAlice) printf("AliRun object found on file\n");
	33	if (!gAlice) {
	34	printf("\n create new gAlice object");
	35	gAlice = new AliRun("gAlice","Alice test program");
	36	}
	37	}
	38
	39	// Create some histograms
	40
	41	TH1F Hits[10], HitDensity[10], Occ0[10],Occ1[10], *Mult[10];
	42	for (Int_t i=0; i<10; i++) {
	43	Hits[i] = new TH1F("hits1","Hit Distribution",26,0,260);
	44	Hits[i]->SetFillColor(0);
	45	Hits[i]->SetXTitle("R (cm)");
	46
	47	HitDensity[i] = new TH1F("dhits1","Hit Density Distribution",26,0,260);
	48	HitDensity[i]->SetFillColor(0);
	49	HitDensity[i]->SetXTitle("R (cm)");
	50
	51	Occ0[i] = new TH1F("occ0","Occupancy Density",26,0,260);
	52	Occ0[i] -> SetFillColor(0);
	53	Occ0[i] -> SetXTitle("R (cm)");
	54	Occ1[i] = new TH1F("occ1","Occupancy",26,0,260);
	55	Occ1[i] -> SetFillColor(0);
	56	Occ1[i] -> SetXTitle("R (cm)");
	57	Occ1[i] -> SetYTitle("Occupancy");
	58	Mult[i] = new TH1F("mult","Mult distribution",26,0,260);
	59	Mult[i] -> SetFillColor(0);
	60	Mult[i] -> SetXTitle("R (cm)");
	61	}
	62
	63
	64	TH1F *theta = new TH1F("theta","Theta distribution",180,0,180);
65
66	TH1F *emult = new TH1F("emult","Event Multiplicity",100,0,1000);
67	AliMUONChamber* iChamber;
68	AliMUONSegmentation* seg;
69
70	//
71	// Loop over events
72	//
73	Float_t dpx[2], dpy[2];
74	Int_t mdig =0;
75	Int_t nhit=0;
76
77	for (Int_t nev=0; nev<= evNumber2; nev++) {
78	Int_t nparticles = gAlice->GetEvent(nev);
79	if (nev < evNumber1) continue;
80	if (nparticles <= 0) return;
81
82	AliMUON MUON = (AliMUON) gAlice->GetModule("MUON");
83	printf("\n track %d %d \n ", nev, MUON);
84
85	TTree *TH = gAlice->TreeH();
86	Int_t ntracks = TH->GetEntries();
87	//
88	// Loop over events
89	//
90	Int_t EvMult=0;
91
92	for (Int_t track=0; track<ntracks;track++) {
93	gAlice->ResetHits();
94	Int_t nbytes += TH->GetEvent(track);
95	if (MUON) {
96	//
97	// Loop over hits
98	//
99	Int_t NperCh=0;
100
101	for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
102	mHit;
103	mHit=(AliMUONHit*)MUON->NextHit())
104	{
bd21d472	105	Int_t nch = mHit->Chamber(); // chamber number
	106	Float_t x = mHit->X(); // x-pos of hit
	107	Float_t y = mHit->Y(); // y-pos
	108	Float_t Eloss = mHit->Eloss();
	109	Float_t Theta = mHit->Theta();
a9e2aefa	110	theta->Fill(Theta,(float) 1);
	111	if (nch >10) continue;
	112	if (nch ==1) EvMult++;
	113	//
	114	//
	115	iChamber = & MUON->Chamber(nch-1);
	116	response= iChamber.GetResponseModel();
	117	seg = iChamber.GetSegmentationModel(1);
	118	NperCh++;
	119	Int_t i,j;
	120	seg->GetPadIxy(x,y,i,j);
	121	Int_t isec = seg->Sector(i,j);
	122	if (isec ==1 && nch==ic) nhit++;
	123
	124	Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
	125	Float_t r=TMath::Sqrt(xx+yy);
	126	Float_t wgt=1/(210TMath::Pi()*r)/(evNumber2+1);
	127	wgt=wgt(1.+24./(2.TMath::Pi()*r));
	128	Hits[nch-1]->Fill(r,(float) 1);
	129	HitDensity[nch-1]->Fill(r,wgt);
	130	Occ0[nch-1]->Fill(r,wgt*a);
	131	} // hit loop
	132	} // if MUON
	133	} // track loop
	134	emult->Fill(Float_t(EvMult), (Float_t) 1);
	135
	136	Int_t iseg=1;
	137
	138	for (Int_t ich=0; ich<10; ich++) {
	139	iChamber = & MUON->Chamber(ich);
	140	seg=iChamber.GetSegmentationModel(iseg);
	141	gAlice->ResetDigits();
	142	gAlice->TreeD()->GetEvent(iseg);
	143	TClonesArray *MUONDigits = MUON->DigitsAddress(ich);
	144	Int_t Ndigits=MUONDigits->GetEntriesFast();
	145	AliMUONDigit* dig;
	146	printf("\n Reading %d digits\n", Ndigits);
	147	if (MUONDigits) {
	148	for (Int_t ndig=0; ndig<Ndigits; ndig++)
	149	{
	150	dig = (AliMUONDigit*)MUONDigits->UncheckedAt(ndig);
bd21d472	151	Int_t i=dig->PadX();
bd21d472	152	Int_t j=dig->PadY();
a9e2aefa	153	Float_t x,y;
	154	seg->GetPadCxy(i,j,x,y);
	155	Int_t isec = seg->Sector(i,j);
	156	Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
	157	Float_t r=TMath::Sqrt(xx+yy);
	158	// if (r<25)
	159	// printf("\n Sector,a %d %f", isec,a);
	160	if (isec==1) mdig++;
	161
	162	Float_t wgt;
	163	wgt=1/(2.10.TMath::Pi()r)/(evNumber2+1)a;
	164	// Take into account inefficiency due to frames
	165	wgt=wgt(1.+24./(2.TMath::Pi()*r));
	166	Occ1[ich]->Fill(r,wgt);
	167	} // digit loop
	168	Mult[ich]->Divide(Occ1[ich],Occ0[ich]);
	169	} // chamber loop
	170	} // if MUONDigits
	171	} // event loop
	172	//
	173	printf("\n hits, digits %d %d\n ", nhit, mdig);
	174
	175
	176
	177	//Create a canvas, set the view range, show histograms
	178	TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
	179	pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
	180	pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
	181	pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
	182	pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
	183	pad11->SetFillColor(0);
	184	pad12->SetFillColor(0);
	185	pad13->SetFillColor(0);
	186	pad14->SetFillColor(0);
	187	pad11->Draw();
	188	pad12->Draw();
	189	pad13->Draw();
	190	pad14->Draw();
	191
	192	pad11->cd();
	193	HitDensity[0]->Draw();
	194	pad12->cd();
	195	HitDensity[1]->Draw();
	196	pad13->cd();
	197	HitDensity[2]->Draw();
	198	pad14->cd();
	199	HitDensity[3]->Draw();
	200
	201	TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
	202	pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
	203	pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
	204	pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
	205	pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
	206	pad21->SetFillColor(0);
	207	pad22->SetFillColor(0);
	208	pad23->SetFillColor(0);
	209	pad24->SetFillColor(0);
	210	pad21->Draw();
	211	pad22->Draw();
	212	pad23->Draw();
	213	pad24->Draw();
	214
	215	pad21->cd();
	216	HitDensity[4]->Draw();
217	pad22->cd();
218	HitDensity[5]->Draw();
219	pad23->cd();
220	HitDensity[6]->Draw();
221	pad24->cd();
222	HitDensity[7]->Draw();
223
224	TCanvas *c3 = new TCanvas("c3","Hit Densities",400,10,600,700);
225	pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
226	pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
227	pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
228	pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
229	pad31->SetFillColor(0);
230	pad32->SetFillColor(0);
231	pad33->SetFillColor(0);
232	pad34->SetFillColor(0);
233	pad31->Draw();
234	pad32->Draw();
235	pad33->Draw();
236	pad34->Draw();
237
238	pad31->cd();
239	HitDensity[8]->Draw();
240	pad32->cd();
241	HitDensity[9]->Draw();
242
243	TCanvas *c4 = new TCanvas("c4","Occupancies",400,10,600,700);
244	pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
245	pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
246	pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
247	pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
248	pad41->SetFillColor(0);
249	pad42->SetFillColor(0);
250	pad43->SetFillColor(0);
251	pad44->SetFillColor(0);
252	pad41->Draw();
253	pad42->Draw();
254	pad43->Draw();
255	pad44->Draw();
256
257
258	pad41->cd();
259	Occ1[0]->Draw();
260	pad42->cd();
261	Occ1[2]->Draw();
262	pad43->cd();
263	Occ1[4]->Draw();
264	pad44->cd();
265	Occ1[6]->Scale(1.25);
266	Occ1[6]->Draw();
267
268	TCanvas *c5 = new TCanvas("c5","Occupancies",400,10,600,700);
269	pad51 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
270	pad52 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
271	pad53 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
272	pad54 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
273	pad51->SetFillColor(0);
274	pad52->SetFillColor(0);
275	pad53->SetFillColor(0);
276	pad54->SetFillColor(0);
277	pad51->Draw();
278	pad52->Draw();
279	pad53->Draw();
280	pad54->Draw();
281
282
283	pad51->cd();
284	Occ1[8]->Scale(1.25);
285
286	Occ1[8]->Draw();
287	}
288
289
290
291
292
293
294
295
296
297