Corrections by G.Chabratova
[u/mrichter/AliRoot.git] / MUON / MUONcombi.C
1 void 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 /*
13 <img src="gif/anal.gif">
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();
40    TParticle *Part;
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;
100    TParticle* Muon1, *Muon2;
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