get tables from the aliroot directory if they are not in the current one
[u/mrichter/AliRoot.git] / MUON / MUONcombi.C
1 void MUONcombi (Int_t evNumber=0)
2 {
3     const Float_t runWeight = 4.e8;
4
5 // Dynamically link some shared libs
6    if (gClassTable->GetID("AliRun") < 0) {
7       gROOT->LoadMacro("loadlibs.C");
8       loadlibs();
9    } else {
10       delete gAlice;
11       gAlice = 0;
12    }
13
14 // Connect the Root Galice file containing Geometry, Kine and Hits
15    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
16    if (!file) file = new TFile("galice.root");
17
18 // Get AliRun object from file or create it if not on file
19    if (!gAlice) {
20       gAlice = (AliRun*)file->Get("gAlice");
21       if (gAlice) printf("AliRun object found on file\n");
22       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
23    }
24
25 // Import the Kine and Hits Trees for the event evNumber in the file
26    Int_t nparticles = gAlice->GetEvent(evNumber);
27    if (nparticles <= 0) return;
28 //
29 //
30    TH1F *dmass  = new TH1F("dmass","Dimuon-Mass Distribution"
31                           ,100,0.,5.);
32    
33    TH1F *dmassc = new TH1F("dmassc","Dimuon-Mass Distribution"
34                            ,200,0.,10.);
35    TH1F *dmassd = new TH1F("dmassd","Dimuon-Mass Distribution"
36                            ,200,0.,10.);
37
38    TH1F *pt     = new TH1F("pt","pT-single"
39                            ,50,0.,10.);
40    TH1F *hCont[30];
41    
42 //
43 //  Generator Loop
44 //
45
46    Int_t i=0;
47    
48    AliGenCocktailEntry *Entry, *e1, *e2;
49    AliGenCocktail*  Cocktail = (AliGenCocktail*) gAlice->Generator();
50 //                   Single Generator
51    for (Entry=Cocktail->FirstGenerator();
52         Entry;
53         Entry=Cocktail->NextGenerator()
54        ) {
55        Entry->PrintInfo();
56      }
57 //
58 // Initialize Combinator 
59 //
60    AliDimuCombinator* Combinator = new AliDimuCombinator();
61    Combinator->SetEtaCut(2.5, 4.);
62    Combinator->SetPtMin(1.0);   
63    
64    i=0;
65    
66 //
67 // Single Muon  Loop  
68 //
69    Combinator->ResetRange();
70    
71    
72
73    for(Muon=Combinator->FirstMuon();
74        Muon; 
75        Muon=Combinator->NextMuon()) {
76 //
77        Int_t chfirst = Muon->GetFirstDaughter();
78        Int_t chlast  = Muon->GetLastDaughter();
79        Int_t parent  = Muon->GetFirstMother();
80        Float_t ptm   = Muon->Pt();
81        Float_t eta   = Muon->Eta();
82        
83 //       printf("\n Particle %d Parent %d first child %d last child %d",
84 //            i,parent, chfirst, chlast);
85 //       printf("\n Particle pt, eta: %f , %f ", ptm, eta);
86        i++;
87    }
88 //
89 // Di-Muon Loop
90
91    Float_t pt1,pt2;
92    TParticle* Muon1;
93    TParticle* Muon2;
94
95    Combinator->ResetRange();
96
97 //
98 //   Dimuon Loop controlled by Generator Loop
99 //
100    Float_t  sig = 0;
101    Int_t  icont = 0;
102    char name[30];   
103    for (Cocktail->FirstGeneratorPair(e1,e2);
104         (e1&&e2);
105         Cocktail->NextGeneratorPair(e1,e2)
106        ){
107
108        sprintf(name, "%s-%s", e1->GetName(), e2->GetName());
109        hCont[icont] = new TH1F(name,"Dimuon-Mass Distribution",100,0.,5.);
110        printf("\n ----------------------------------------------------");
111        e1->PrintInfo();
112        e2->PrintInfo();
113        printf("\n ----------------------------------------------------");
114        Combinator->SetFirstRange (e1->GetFirst(), e1->GetLast());
115        Combinator->SetSecondRange(e2->GetFirst(), e2->GetLast());
116        Combinator->SetRate(e1->Rate(), e2->Rate());
117        for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
118             (Muon1 && Muon2);
119             Combinator->NextMuonPairSelected(Muon1,Muon2))
120        {
121            pt1 = Muon1->Pt();
122            pt2 = Muon2->Pt();       
123            Float_t mass = Combinator->Mass(Muon1, Muon2);
124            Float_t wgt  = runWeight*Combinator->Weight(Muon1, Muon2)*
125                Combinator->DecayProbability(Muon1)*
126                Combinator->DecayProbability(Muon2);
127            pt->Fill(pt1, wgt);
128            pt->Fill(pt2, wgt);
129            Float_t smeared_mass = mass;
130            Combinator->SmearGauss(0.02*mass, smeared_mass);
131            
132            if (TMath::Min(pt1,pt2) > -0.5*TMath::Max(pt1,pt2)+2.) {
133                if (Combinator->Correlated(Muon1, Muon2)) {
134                    dmassc->Fill(smeared_mass, wgt);
135                    sig += wgt;
136                } else {
137 // account for the fact that we sum like-sign and unlike-sign              
138                    wgt *= 0.5;
139                }
140                dmass->Fill(smeared_mass, wgt);
141                hCont[icont]->Fill(smeared_mass, wgt);
142            }
143        } // Dimuon Loop
144        icont++;
145    }// Generator Loop
146    printf("\n Signal %e \n \n", sig);
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    pad12->SetLogy(0);
172    dmass->SetFillColor(42);
173    dmass->SetXTitle("m (GeV)");
174    dmass->Draw();
175
176    pad13->cd();
177    pad13->SetLogy(0);
178    dmassc->SetFillColor(42);
179    dmassc->SetXTitle("m (GeV)");
180    dmassc->Draw();
181
182    pad14->cd();
183    pad14->SetLogy(0);
184    dmassd->SetFillColor(42);
185    dmassd->SetXTitle("m (GeV)");
186    dmassd->Draw();
187
188
189 }
190
191
192
193
194
195
196
197
198
199
200