Corrections by G.Chabratova
[u/mrichter/AliRoot.git] / MUON / MUONcombi.C
CommitLineData
fe4da5cc 1void MUONcombi (Int_t evNumber=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 fill some histograms.
7//
8// Root > .L anal.C //this loads the macro in memory
9// Root > anal(); //by default process first event
10// Root > anal(2); //process third event
11//Begin_Html
12/*
a897a37a 13<img src="gif/anal.gif">
fe4da5cc 14*/
15//End_Html
16/////////////////////////////////////////////////////////////////////////
17// Dynamically link some shared libs
18 if (gClassTable->GetID("AliRun") < 0) {
19 gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
20 gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
21 gSystem->Load("$ROOTSYS/lib/libEG.so");
22 gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
23 gSystem->Load("libGeant3Dummy.so"); //a dummy version of Geant3
24 gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes
25 gSystem->Load("libgalice.so"); // the standard Alice classes
26 }
27//
28// Connect the Root Galice file containing Geometry, Kine and Hits
29 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
30 if (!file) file = new TFile("galice.root");
31
32// Get AliRun object from file or create it if not on file
33 if (!gAlice) {
34 gAlice = (AliRun*)file->Get("gAlice");
35 if (gAlice) printf("AliRun object found on file\n");
36 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
37 }
38// TClonesArray &Partarray = *(gAlice->Particles());
39 TClonesArray *PartArray = gAlice->Particles();
50f986db 40 TParticle *Part;
fe4da5cc 41
42
43// Import the Kine and Hits Trees for the event evNumber in the file
44 Int_t nparticles = gAlice->GetEvent(evNumber);
45 if (nparticles <= 0) return;
46//
47//
48 TH1F *dmass = new TH1F("dmass","Dimuon-Mass Distribution"
49 ,50,0.,10.);
50 TH1F *pt = new TH1F("pt","pT-single"
51 ,50,0.,10.);
52//
53// Generator Loop
54//
55
56
57 AliGenCocktailEntry *Entry, *e1, *e2;
58 AliGenCocktail* Cocktail = (AliGenCocktail*) gAlice->Generator();
59// Single Generator
60 for (Entry=Cocktail->FirstGenerator();
61 Entry;
62 Entry=Cocktail->NextGenerator()
63 ) {
64 Entry->PrintInfo();
65 }
66// Pairs of Generators
67
68//
69// Initialize Combinator
70 DimuonCombinator Combinator = DimuonCombinator(PartArray);
71 Combinator->SetEtaCut(-4, 4);
72 Combinator->SetPtMin(0.5);
73
74 Int_t i=0;
75//
76// Single Muon Loop
77//
78/*
79
80 for(Muon=Combinator.FirstMuon();
81 Muon;
82 Muon=Combinator.NextMuon()) {
83//
84 Int_t chfirst= Muon->GetFirstChild();
85 Int_t chlast = Muon->GetLastChild();
86 Int_t parent = Muon->GetParent();
87 Float_t ptm = Muon->GetPT();
88 Float_t eta = Muon->GetEta();
89
90 printf("\n Particle %d Parent %d first child %d last child %d",
91 i,parent, chfirst, chlast);
92 printf("\n Particle pt, eta: %f , %f ",pt,eta);
93 i++;
94 }
95*/
96//
97// Di-Muon Loop
98
99 Float_t pt1,pt2;
50f986db 100 TParticle* Muon1, *Muon2;
fe4da5cc 101/*
102 Combinator->ResetRange();
103 for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
104 (Muon1 && Muon2);
105 Combinator->NextMuonPairSelected(Muon1,Muon2))
106 {
107 pt1=Muon1->GetPT();
108 pt2=Muon2->GetPT();
109 Float_t mass=Combinator->Mass(Muon1, Muon2);
110 Float_t wgt =Combinator->Weight(Muon1, Muon2);
111 pt->Fill(pt1, wgt);
112 pt->Fill(pt2, wgt);
113 Float_t smeared_mass=mass;
114 Combinator->SmearGauss(0.05*mass, smeared_mass);
115 dmass->Fill(smeared_mass , wgt);
116 }
117*/
118//
119// Dimuon Loop controlled by Generator Loop
120//
121 for (Cocktail->FirstGeneratorPair(e1,e2);
122 (e1&&e2);
123 Cocktail->NextGeneratorPair(e1,e2)
124 ){
125 printf("\n ----------------------------------------------------");
126 e1->PrintInfo();
127 e2->PrintInfo();
128 printf("\n ----------------------------------------------------");
129 Combinator->SetFirstRange (e1->GetFirst(), e1->GetLast());
130 Combinator->SetSecondRange(e2->GetFirst(), e2->GetLast());
131 Combinator->SetRate(e1->Rate(), e2->Rate());
132 for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
133 (Muon1 && Muon2);
134 Combinator->NextMuonPairSelected(Muon1,Muon2))
135 {
136 pt1=Muon1->GetPT();
137 pt2=Muon2->GetPT();
138 Float_t mass=Combinator->Mass(Muon1, Muon2);
139 Float_t wgt =Combinator->Weight(Muon1, Muon2);
140 pt->Fill(pt1, wgt);
141 pt->Fill(pt2, wgt);
142 Float_t smeared_mass=mass;
143 Combinator->SmearGauss(0.05*mass, smeared_mass);
144 dmass->Fill(smeared_mass , wgt);
145 } // Dimuon Loop
146 }// Generator Loop
147
148//
149//Create a canvas, set the view range, show histograms
150//
151 TCanvas *c1 = new TCanvas("c1","Dimuon Plots",400,10,600,700);
152 TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
153 TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
154 TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
155 TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
156 pad11->SetFillColor(11);
157 pad12->SetFillColor(11);
158 pad13->SetFillColor(11);
159 pad14->SetFillColor(11);
160 pad11->Draw();
161 pad12->Draw();
162 pad13->Draw();
163 pad14->Draw();
164
165 pad11->cd();
166 pt->SetFillColor(42);
167 pt->SetXTitle("pT (GeV)");
168 pt->Draw();
169
170 pad12->cd();
171 dmass->SetFillColor(42);
172 dmass->SetXTitle("m (GeV)");
173 dmass->Draw();
174}
175
176
177
178
179
180
181
182
183
184