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

void MUONtestabso (Int_t evNumber1=0,Int_t evNumber2=0) 
{
/////////////////////////////////////////////////////////////////////////
//   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                    
    static Float_t xmuon, ymuon;
    
   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 *zv = new TH1F("zv","z vertex" ,140,470.,505.);
   TH1F *xv = new TH1F("xv","x vertex" ,140,-70.,70.);
   TH1F *yv = new TH1F("yv","y vertex", 140,-70.,70.);

   TH1F *ip  = new TH1F("ip","geant part",50,0.,10.);
   TH2F *rzv  = new TH2F("rzv","R-z vert",100,502,504.,100,0.,100.);

   TH1F *ptUps  = new TH1F("ptUps","pT Upsilon",50,0.,25.);
   TH1F *hde  = new TH1F("hde","dE",200,0.,50);
   TH1F *hde2  = new TH1F("hde2","dE",100,1.5,5.5);
   TH2F *hdevsn  = new TH2F("hdevsn","dE vs N electron",100,0.,15., 20, 0.,20.);

   TH1F *ekine  = new TH1F("ekine","E_kin electrons",70,-5,2);
   TH1F *etheta = new TH1F("etheta","Theta electrons",90,0,90);
   TH1F *edr = new TH1F("edr","Distance to muon",100,0,10);

   TH1F *de  = new TH1F("de","correction",100,-1,1);

   TH1F *dtheta = new TH1F("dtheta","Delta Theta" ,200,-5.,5.);
   AliMUONChamber*  iChamber;
//
//   Loop over events 
//
   Int_t Nh=0;
   Int_t Nh1=0;
       Int_t Nel1=0;
       Int_t Nel2=0;
       Int_t Nel3=0;
       Int_t Nel4=0;
       
   for (Int_t nev=0; nev<= evNumber2; nev++) {
       //cout << "nev         " << nev <<endl;
       Int_t nparticles = gAlice->GetEvent(nev);
       //cout << "nparticles  " << nparticles <<endl;
       if (nev < evNumber1) continue;
       if (nparticles <= 0) return;
       
       AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");


       TTree *TH = gAlice->TreeH();
       Int_t ntracks = TH->GetEntries();
//
//   Loop over tracks
//

       Float_t dE;
       
       for (Int_t track=0; track<ntracks;track++) {
	   Int_t Nelt=0;	   
	   printf("\n nev %d %d %d %d %d\n ", nev, Nel1, Nel2, Nel3, Nel4);
	   gAlice->ResetHits();
	   Int_t nbytes += TH->GetEvent(track);
	   if (MUON)  {
	       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 z     = mHit->fZ;        // y-pos
		   
		   if (nch != 1) continue;
		   
		   Int_t ipart = mHit->fParticle;
		   TClonesArray *fPartArray = gAlice->Particles();
		   TParticle *Part;
		   Int_t ftrack = mHit->fTrack;
		   Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
		   Int_t ipart = Part->GetPdgCode();
		   ip->Fill((float)ipart);
		   TParticle *Mother;
		   
		   Float_t px0=Part->Px();
		   Float_t py0=Part->Py();		   
		   Float_t pz0=Part->Pz();
		   Float_t thetax0=TMath::ATan2(px0,pz0);
		   Float_t thetay0=TMath::ATan2(py0,pz0);		   
		   
		   if (ipart == kMuonPlus || ipart == kMuonMinus) {
//		       Int_t imo = Part->GetFirstMother();
//		       Mother = (TParticle*) fPartArray->UncheckedAt(imo);
//		       Float_t pt = Mother->Pt();
//		       ptUps->Fill(pt, (float) 1);
		       Float_t E=Part->Energy();
		       Float_t Eloc=mHit->fPTot;
		       Float_t corr=Eloc+CorrectP(Eloc,mHit->fTheta);
		       printf("\n %f %f %f", E, Eloc, corr);
		       de->Fill(E-corr,1.);
		       dE = E-Eloc;
		       if (dE<50)  hde->Fill(dE, (float) 1);
		       if (dE<5.5) hde2->Fill(dE, (float) 1);
		       xmuon=mHit->fX;
		       ymuon=mHit->fY;
		       Float_t thetax=TMath::ATan2(mHit->fCxHit, mHit->fPTot);
		       Float_t thetay=TMath::ATan2(mHit->fCyHit, mHit->fPTot);
		       dtheta->Fill((thetax-thetax0)*1000., 1.);
		       dtheta->Fill((thetay-thetay0)*1000., 1.);		       
		   }
		   
		   if (ipart == kElectron || ipart == kPositron) {

		   
		       Float_t xvert = Part->Vx();       //  vertex  
		       Float_t yvert = Part->Vy();       //  vertex
		       Float_t zvert = Part->Vz();       // z vertex 
		       if (zvert < 503 && mHit->fTheta<90) {
			   Nelt++;
			   Float_t px = Part->Px();
			   Float_t py = Part->Py();
			   Float_t pz = Part->Pz();
			   Float_t Ek = Part->Energy()-Part->GetMass();
			   
			   Int_t imo = Part->GetFirstMother();
			   Mother = (TParticle*) fPartArray->UncheckedAt(imo);
			   Int_t imot = Mother->GetPdgCode();
			   xv->Fill(xvert);
			   yv->Fill(yvert);
			   zv->Fill(zvert);
			   Float_t rvert=TMath::Sqrt(xvert*xvert+yvert*yvert);
			   rzv->Fill(zvert,rvert);
			   ekine->Fill(TMath::Log10(Ek),1.);
			   etheta->Fill(mHit->fTheta,1.);
			   Float_t ex=mHit->fX;
			   Float_t ey=mHit->fY;
			   dr=TMath::Sqrt((ex-xmuon)*(ex-xmuon)+(ey-ymuon)*(ey-ymuon));
			   edr->Fill(dr,1.);
			   
		       }
		   }
		   
	       } // hits 
	   } // if MUON 
	   if (Nelt == 1) Nel1++;
	   if (Nelt == 2) Nel2++;
	   if (Nelt == 3) Nel3++;
	   if (Nelt  > 3) Nel4++;
	   hdevsn->Fill(dE, (float) Nelt, (float) 1);

       } // tracks

    } // event

   
//Create a canvas, set the view range, show histograms
   TCanvas *c1 = new TCanvas("c1","Vetices from electrons and positrons",400,10,600,700);
   c1->Divide(2,2);
   
   c1->cd(1);
   ip->SetFillColor(42);
   ip->SetXTitle("ipart");
   ip->Draw();

   c1->cd(2);
   xv->SetFillColor(42);
   xv->SetXTitle("xvert");
   xv->Draw();

   c1->cd(3);
   yv->SetFillColor(42);
   yv->SetXTitle("yvert");
   yv->Draw();

   c1->cd(4);
   zv->SetFillColor(42);
   zv->SetXTitle("zvert");
   zv->Draw();

   TCanvas *c2 = new TCanvas("c2"," ",400,10,600,700);
   c2->Divide(2,2);
   c2->cd(1);
   rzv->SetXTitle("zvert");
   rzv->SetYTitle("rvert");
   rzv->Draw();

   c2->cd(2);
   ptUps->SetXTitle("pt");
   ptUps->Draw();

   c2->cd(3);
   hde->SetXTitle("dE");
   hde->Draw();


   c2->cd(4);
   hde2->SetXTitle("dE");
   hde2->Draw();
   TCanvas *c3 = new TCanvas("c3"," ",400,10,600,700);
   c3->Divide(2,2);
   c3->cd(1);
   ekine->SetXTitle("E_kin");
   ekine->Draw();

   c3->cd(2);
   etheta->SetXTitle("Theta");
   etheta->Draw();

   c3->cd(3);
   edr->SetXTitle("Distance to muon");
   edr->Draw();

   c3->cd(4);
   dtheta->SetXTitle(" ");
   dtheta->Draw();

}

Float_t CorrectP(Float_t p, Float_t theta)
{
    printf("\n %f%", theta);
    
    if (theta<3.) {
//W
	if (p<15) {
	    return 2.737+0.0494*p-0.001123*p*p;
	} else {
	    return 3.0643+0.01346*p;
	}
    } else {
//Pb
	if (p<15) {
	    return 2.1380+0.0351*p-0.000853*p*p;
	} else {
	    return 2.407+0.00702*p;
	}
    }
}
Commit	Line	Data
a897a37a	1	void MUONtestabso (Int_t evNumber1=0,Int_t evNumber2=0)
	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
a9e2aefa	11	static Float_t xmuon, ymuon;
a9e2aefa	12
a897a37a	13	if (gClassTable->GetID("AliRun") < 0) {
	14	gROOT->LoadMacro("loadlibs.C");
	15	loadlibs();
	16	}
	17
	18	// Connect the Root Galice file containing Geometry, Kine and Hits
	19
	20	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject("galice.root");
	21
	22
	23	if (!file) {
	24	printf("\n Creating galice.root \n");
	25	file = new TFile("galice.root");
	26	} else {
	27	printf("\n galice.root found in file list");
	28	}
	29	file->ls();
	30
	31	// Get AliRun object from file or create it if not on file
	32	if (!gAlice) {
	33	gAlice = (AliRun*)(file->Get("gAlice"));
	34	if (gAlice) printf("AliRun object found on file\n");
	35	if (!gAlice) {
a9e2aefa	36	printf("\n Create new gAlice object");
a897a37a	37	gAlice = new AliRun("gAlice","Alice test program");
	38	}
	39	}
	40
a897a37a	41
	42	// Create some histograms
	43
	44
a9e2aefa	45	TH1F *zv = new TH1F("zv","z vertex" ,140,470.,505.);
	46	TH1F *xv = new TH1F("xv","x vertex" ,140,-70.,70.);
	47	TH1F *yv = new TH1F("yv","y vertex", 140,-70.,70.);
	48
	49	TH1F *ip = new TH1F("ip","geant part",50,0.,10.);
	50	TH2F *rzv = new TH2F("rzv","R-z vert",100,502,504.,100,0.,100.);
	51
	52	TH1F *ptUps = new TH1F("ptUps","pT Upsilon",50,0.,25.);
	53	TH1F *hde = new TH1F("hde","dE",200,0.,50);
	54	TH1F *hde2 = new TH1F("hde2","dE",100,1.5,5.5);
	55	TH2F *hdevsn = new TH2F("hdevsn","dE vs N electron",100,0.,15., 20, 0.,20.);
a897a37a	56
a9e2aefa	57	TH1F *ekine = new TH1F("ekine","E_kin electrons",70,-5,2);
	58	TH1F *etheta = new TH1F("etheta","Theta electrons",90,0,90);
	59	TH1F *edr = new TH1F("edr","Distance to muon",100,0,10);
a897a37a	60
a9e2aefa	61	TH1F *de = new TH1F("de","correction",100,-1,1);
	62
	63	TH1F *dtheta = new TH1F("dtheta","Delta Theta" ,200,-5.,5.);
	64	AliMUONChamber* iChamber;
a897a37a	65	//
	66	// Loop over events
	67	//
	68	Int_t Nh=0;
	69	Int_t Nh1=0;
a9e2aefa	70	Int_t Nel1=0;
	71	Int_t Nel2=0;
	72	Int_t Nel3=0;
	73	Int_t Nel4=0;
	74
a897a37a	75	for (Int_t nev=0; nev<= evNumber2; nev++) {
a9e2aefa	76	//cout << "nev " << nev <<endl;
a897a37a	77	Int_t nparticles = gAlice->GetEvent(nev);
a9e2aefa	78	//cout << "nparticles " << nparticles <<endl;
a897a37a	79	if (nev < evNumber1) continue;
a897a37a	80	if (nparticles <= 0) return;
a9e2aefa	81
a897a37a	82	AliMUON MUON = (AliMUON) gAlice->GetModule("MUON");
a9e2aefa	83
a897a37a	84
	85	TTree *TH = gAlice->TreeH();
	86	Int_t ntracks = TH->GetEntries();
	87	//
	88	// Loop over tracks
	89	//
a9e2aefa	90
	91	Float_t dE;
	92
a897a37a	93	for (Int_t track=0; track<ntracks;track++) {
a9e2aefa	94	Int_t Nelt=0;
a9e2aefa	95	printf("\n nev %d %d %d %d %d\n ", nev, Nel1, Nel2, Nel3, Nel4);
a897a37a	96	gAlice->ResetHits();
	97	Int_t nbytes += TH->GetEvent(track);
	98	if (MUON) {
a9e2aefa	99	for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
a897a37a	100	mHit;
a9e2aefa	101	mHit=(AliMUONHit*)MUON->NextHit())
a897a37a	102	{
	103	Int_t nch = mHit->fChamber; // chamber number
	104	Float_t x = mHit->fX; // x-pos of hit
	105	Float_t y = mHit->fY; // y-pos
	106	Float_t z = mHit->fZ; // y-pos
a9e2aefa	107
	108	if (nch != 1) continue;
	109
	110	Int_t ipart = mHit->fParticle;
	111	TClonesArray *fPartArray = gAlice->Particles();
	112	TParticle *Part;
	113	Int_t ftrack = mHit->fTrack;
	114	Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
	115	Int_t ipart = Part->GetPdgCode();
	116	ip->Fill((float)ipart);
	117	TParticle *Mother;
	118
	119	Float_t px0=Part->Px();
	120	Float_t py0=Part->Py();
	121	Float_t pz0=Part->Pz();
	122	Float_t thetax0=TMath::ATan2(px0,pz0);
	123	Float_t thetay0=TMath::ATan2(py0,pz0);
	124
	125	if (ipart == kMuonPlus \|\| ipart == kMuonMinus) {
	126	// Int_t imo = Part->GetFirstMother();
	127	// Mother = (TParticle*) fPartArray->UncheckedAt(imo);
	128	// Float_t pt = Mother->Pt();
	129	// ptUps->Fill(pt, (float) 1);
	130	Float_t E=Part->Energy();
	131	Float_t Eloc=mHit->fPTot;
	132	Float_t corr=Eloc+CorrectP(Eloc,mHit->fTheta);
	133	printf("\n %f %f %f", E, Eloc, corr);
	134	de->Fill(E-corr,1.);
	135	dE = E-Eloc;
	136	if (dE<50) hde->Fill(dE, (float) 1);
	137	if (dE<5.5) hde2->Fill(dE, (float) 1);
	138	xmuon=mHit->fX;
	139	ymuon=mHit->fY;
	140	Float_t thetax=TMath::ATan2(mHit->fCxHit, mHit->fPTot);
	141	Float_t thetay=TMath::ATan2(mHit->fCyHit, mHit->fPTot);
	142	dtheta->Fill((thetax-thetax0)*1000., 1.);
	143	dtheta->Fill((thetay-thetay0)*1000., 1.);
	144	}
	145
	146	if (ipart == kElectron \|\| ipart == kPositron) {
	147
	148
	149	Float_t xvert = Part->Vx(); // vertex
	150	Float_t yvert = Part->Vy(); // vertex
	151	Float_t zvert = Part->Vz(); // z vertex
	152	if (zvert < 503 && mHit->fTheta<90) {
	153	Nelt++;
	154	Float_t px = Part->Px();
	155	Float_t py = Part->Py();
	156	Float_t pz = Part->Pz();
	157	Float_t Ek = Part->Energy()-Part->GetMass();
	158
	159	Int_t imo = Part->GetFirstMother();
	160	Mother = (TParticle*) fPartArray->UncheckedAt(imo);
	161	Int_t imot = Mother->GetPdgCode();
	162	xv->Fill(xvert);
	163	yv->Fill(yvert);
	164	zv->Fill(zvert);
	165	Float_t rvert=TMath::Sqrt(xvertxvert+yvertyvert);
	166	rzv->Fill(zvert,rvert);
	167	ekine->Fill(TMath::Log10(Ek),1.);
	168	etheta->Fill(mHit->fTheta,1.);
	169	Float_t ex=mHit->fX;
	170	Float_t ey=mHit->fY;
171	dr=TMath::Sqrt((ex-xmuon)(ex-xmuon)+(ey-ymuon)(ey-ymuon));
172	edr->Fill(dr,1.);
173
174	}
175	}
176
177	} // hits
178	} // if MUON
179	if (Nelt == 1) Nel1++;
180	if (Nelt == 2) Nel2++;
181	if (Nelt == 3) Nel3++;
182	if (Nelt > 3) Nel4++;
183	hdevsn->Fill(dE, (float) Nelt, (float) 1);
184
185	} // tracks
186
187	} // event
a897a37a	188
	189
	190	//Create a canvas, set the view range, show histograms
	191	TCanvas *c1 = new TCanvas("c1","Vetices from electrons and positrons",400,10,600,700);
a9e2aefa	192	c1->Divide(2,2);
a897a37a	193
a9e2aefa	194	c1->cd(1);
a897a37a	195	ip->SetFillColor(42);
	196	ip->SetXTitle("ipart");
	197	ip->Draw();
	198
a9e2aefa	199	c1->cd(2);
a897a37a	200	xv->SetFillColor(42);
	201	xv->SetXTitle("xvert");
	202	xv->Draw();
	203
a9e2aefa	204	c1->cd(3);
a897a37a	205	yv->SetFillColor(42);
	206	yv->SetXTitle("yvert");
	207	yv->Draw();
	208
a9e2aefa	209	c1->cd(4);
a897a37a	210	zv->SetFillColor(42);
	211	zv->SetXTitle("zvert");
	212	zv->Draw();
	213
a9e2aefa	214	TCanvas *c2 = new TCanvas("c2"," ",400,10,600,700);
	215	c2->Divide(2,2);
	216	c2->cd(1);
a897a37a	217	rzv->SetXTitle("zvert");
	218	rzv->SetYTitle("rvert");
	219	rzv->Draw();
a9e2aefa	220
	221	c2->cd(2);
	222	ptUps->SetXTitle("pt");
	223	ptUps->Draw();
	224
	225	c2->cd(3);
	226	hde->SetXTitle("dE");
	227	hde->Draw();
	228
	229
	230	c2->cd(4);
	231	hde2->SetXTitle("dE");
	232	hde2->Draw();
	233	TCanvas *c3 = new TCanvas("c3"," ",400,10,600,700);
	234	c3->Divide(2,2);
	235	c3->cd(1);
	236	ekine->SetXTitle("E_kin");
	237	ekine->Draw();
	238
	239	c3->cd(2);
	240	etheta->SetXTitle("Theta");
	241	etheta->Draw();
	242
	243	c3->cd(3);
	244	edr->SetXTitle("Distance to muon");
	245	edr->Draw();
	246
	247	c3->cd(4);
	248	dtheta->SetXTitle(" ");
	249	dtheta->Draw();
	250
	251	}
	252
	253	Float_t CorrectP(Float_t p, Float_t theta)
	254	{
	255	printf("\n %f%", theta);
	256
	257	if (theta<3.) {
	258	//W
	259	if (p<15) {
	260	return 2.737+0.0494p-0.001123p*p;
	261	} else {
	262	return 3.0643+0.01346*p;
	263	}
	264	} else {
	265	//Pb
	266	if (p<15) {
	267	return 2.1380+0.0351p-0.000853p*p;
	268	} else {
	269	return 2.407+0.00702*p;
	270	}
	271	}
a897a37a	272	}
a9e2aefa	273