/////////////////////////////////////////////////////////////////////////
// Dynamically link some shared libs
-
+ static Float_t xmuon, ymuon;
+
if (gClassTable->GetID("AliRun") < 0) {
gROOT->LoadMacro("loadlibs.C");
loadlibs();
gAlice = (AliRun*)(file->Get("gAlice"));
if (gAlice) printf("AliRun object found on file\n");
if (!gAlice) {
- printf("\n create new gAlice object");
+ printf("\n Create new gAlice object");
gAlice = new AliRun("gAlice","Alice test program");
}
}
- printf ("I'm after gAlice \n");
// Create some histograms
- TH1F *zv = new TH1F("zv","z vertex"
- ,100,450..,506.);
- TH1F *xv = new TH1F("xv","x vertex"
- ,140,-70.,70.);
- TH1F *yv = new TH1F("yv","y vertex"
- ,140,-70.,70.);
+ 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 *ip = new TH1F("ip","geant part"
- ,50,0.,10.);
- TH2F *rzv = new TH2F("rzv","R-z vert"
- ,100,500.,506.,50,0.,50.);
+ 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);
- AliMUONchamber* iChamber;
+ 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;
+ //cout << "nev " << nev <<endl;
Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nparticles " << nparticles <<endl;
+ //cout << "nparticles " << nparticles <<endl;
if (nev < evNumber1) continue;
if (nparticles <= 0) return;
-
+
AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- printf("\n nev %d \n ", nev);
+
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);
+ for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
mHit;
- mHit=(AliMUONhit*)MUON->NextHit())
+ 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 id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
- ip->Fill((float)ipart);
-
- if (ipart > 3) continue;
-
- Float_t xvert = Part->Vx(); // vertex
- Float_t yvert = Part->Vy(); // vertex
- Float_t zvert = Part->Vz(); // z vertex
- xv->Fill(xvert);
- yv->Fill(yvert);
- zv->Fill(zvert);
- Float_t rvert=TMath::Sqrt(xvert*xvert+yvert*yvert);
- rzv->Fill(zvert,rvert);
-
- }
-
- }
- }
- }
+
+ 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->Px(), mHit->Momentum());
+ Float_t thetay=TMath::ATan2(mHit->Py(), mHit->Momentum());
+ 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);
- 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(11);
- pad12->SetFillColor(11);
- pad13->SetFillColor(11);
- pad14->SetFillColor(11);
- pad11->Draw();
- pad12->Draw();
- pad13->Draw();
- pad14->Draw();
+ c1->Divide(2,2);
- pad11->cd();
+ c1->cd(1);
ip->SetFillColor(42);
ip->SetXTitle("ipart");
ip->Draw();
- pad12->cd();
+ c1->cd(2);
xv->SetFillColor(42);
xv->SetXTitle("xvert");
xv->Draw();
- pad13->cd();
+ c1->cd(3);
yv->SetFillColor(42);
yv->SetXTitle("yvert");
yv->Draw();
- pad14->cd();
+ c1->cd(4);
zv->SetFillColor(42);
zv->SetXTitle("zvert");
zv->Draw();
- TCanvas *c2 = new TCanvas("c2","R-Z vertex distribution",400,10,600,700);
+ 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;
+ }
+ }
}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff