b6ee680f0c816591520d4aa090878c0edb62d89d
[u/mrichter/AliRoot.git] / MUON / MUONstraggling.C
1    TH1F *mass1 = new TH1F("mass1","Invariant Mass",120,0,12);
2    TH1F *mass2 = new TH1F("mass2","Invariant Mass",100,9.0,10.5);
3    TH1F *mass3 = new TH1F("mass3","Invariant Mass",100,2.5,3.5);
4 void MUONstraggling (Int_t evNumber1=0,Int_t evNumber2=0) 
5 {
6 /////////////////////////////////////////////////////////////////////////
7 //   This macro is a small example of a ROOT macro
8 //   illustrating how to read the output of GALICE
9 //   and do some analysis.
10 //   
11 /////////////////////////////////////////////////////////////////////////
12
13 // Dynamically link some shared libs                    
14     static Float_t xmuon, ymuon;
15     
16    if (gClassTable->GetID("AliRun") < 0) {
17       gROOT->LoadMacro("loadlibs.C");
18       loadlibs();
19    }
20
21 // Connect the Root Galice file containing Geometry, Kine and Hits
22
23     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
24
25     
26     if (!file) {
27         printf("\n Creating galice.root \n");
28         file = new TFile("galice.root");
29     } else {
30         printf("\n galice.root found in file list");
31     }
32     file->ls();
33     
34 // Get AliRun object from file or create it if not on file
35     if (!gAlice) {
36         gAlice = (AliRun*)(file->Get("gAlice"));
37         if (gAlice) printf("AliRun object found on file\n");
38         if (!gAlice) {
39             printf("\n Create new gAlice object");
40             gAlice = new AliRun("gAlice","Alice test program");
41         }
42     }
43
44     
45 //  Create some histograms
46
47
48    AliMUONChamber*  iChamber;
49 //
50 //   Loop over events 
51 //
52    Int_t Nh=0;
53    Int_t Nh1=0;
54    for (Int_t nev=0; nev<= evNumber2; nev++) {
55        Int_t nparticles = gAlice->GetEvent(nev);
56        //cout << "nparticles  " << nparticles <<endl;
57        if (nev < evNumber1) continue;
58        if (nparticles <= 0) return;
59        
60        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
61
62
63        TTree *TH = gAlice->TreeH();
64        Int_t ntracks = TH->GetEntries();
65 //
66 //   Loop over tracks
67 //
68
69        Float_t pups[4]={0,0,0,0};
70        for (Int_t track=0; track<ntracks;track++) {
71            Int_t Nelt=0;           
72            gAlice->ResetHits();
73            Int_t nbytes += TH->GetEvent(track);
74
75            
76            if (MUON)  {
77                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
78                    mHit;
79                    mHit=(AliMUONHit*)MUON->NextHit()) 
80                {
81                    Int_t   nch   = mHit->Chamber();  // chamber number
82                    Float_t x     = mHit->X();        // x-pos of hit
83                    Float_t y     = mHit->Y();        // y-pos
84                    Float_t z     = mHit->Z();        // y-pos
85                    Float_t p=mHit->Momentum();
86                    Float_t px=mHit->Px();
87                    Float_t py=mHit->Py();
88                    Float_t pz=mHit->Pz();
89                    
90                    if (nch != 1) continue;
91                    
92                    Int_t ipart = mHit->Particle();
93                    TParticle *Part;
94                    Int_t ftrack = mHit->Track();
95                    Part = gAlice->Particle(ftrack);
96                    Int_t ipart = Part->GetPdgCode();
97                    TParticle *Mother;
98                    Float_t px0=Part->Px();
99                    Float_t py0=Part->Py();                 
100                    Float_t pz0=Part->Pz();
101                    Float_t e0=Part->Energy();
102                    
103                    if (ipart == kMuonPlus || ipart == kMuonMinus) {
104 //
105 // Branson Correction
106 //
107                        Float_t zch=505.;
108                        Float_t r=TMath::Sqrt(x*x+y*y);
109                        Float_t zb;
110                        
111                        if (r<26.3611) {
112                            zb=466.;
113                        } else {
114                            zb=441.;
115                        }
116                        
117                        Float_t xb=x-(zch-zb)*px/pz;
118                        Float_t yb=y-(zch-zb)*py/pz;                    
119                        Float_t tx=xb/zb;
120                        Float_t ty=yb/zb;
121                        pz=zb*p/TMath::Sqrt(zb*zb+xb*xb+yb*yb);
122                        px=pz*tx;
123                        py=pz*ty;
124                        
125 //
126 // Energy Correction
127 //
128 //
129                        Float_t corr=(p+CorrectP(p,mHit->fTheta))/p;
130                        pups[0]+=p*corr;
131                        pups[1]+=px*corr;
132                        pups[2]+=py*corr;
133                        pups[3]+=pz*corr;
134                    }
135                } // hits 
136            } // if MUON 
137        } // tracks
138        Float_t mass=TMath::Sqrt(pups[0]*pups[0]-pups[1]*pups[1]-pups[2]*pups[2]-pups[3]*pups[3]);
139        mass1->Fill(mass, 1.);
140        mass2->Fill(mass, 1.);
141        mass3->Fill(mass, 1.);          
142     } // event
143 //Create a canvas, set the view range, show histograms
144    TCanvas *c1 = new TCanvas("c1","Vertices from electrons and positrons",400,10,600,700);
145    c1->Divide(2,2);
146    
147    c1->cd(1);
148    mass1->SetFillColor(42);
149    mass1->SetXTitle("Mass (GeV)");
150    mass1->Draw();
151
152    c1->cd(2);
153    mass2->SetFillColor(42);
154    mass2->SetXTitle("Mass (GeV)");
155    mass2->Draw();
156
157    c1->cd(3);
158    mass3->SetFillColor(42);
159    mass3->SetXTitle("Mass (GeV)");
160    mass3->Draw();
161
162 }
163
164
165
166 Float_t CorrectP(Float_t p, Float_t theta)
167 {
168     if (theta<3.) {
169 //W
170         if (p<15) {
171             return 2.737+0.0494*p-0.001123*p*p;
172         } else {
173             return 3.0643+0.01346*p;
174         }
175     } else {
176 //Pb
177         if (p<15) {
178             return 2.1380+0.0351*p-0.000853*p*p;
179         } else {
180             return 2.407+0.00702*p;
181         }
182     }
183 }
184
185
186
187
188
189