Energy correction for default absorber configuration.
[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 // 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) {
37             printf("\n Create new gAlice object");
38             gAlice = new AliRun("gAlice","Alice test program");
39         }
40     }
41
42     
43 //  Create some histograms
44
45
46    AliMUONChamber*  iChamber;
47 //
48 //   Loop over events 
49 //
50    Int_t Nh=0;
51    Int_t Nh1=0;
52    for (Int_t nev=0; nev<= evNumber2; nev++) {
53        Int_t nparticles = gAlice->GetEvent(nev);
54        //cout << "nparticles  " << nparticles <<endl;
55        if (nev < evNumber1) continue;
56        if (nparticles <= 0) return;
57        
58        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
59
60
61        TTree *TH = gAlice->TreeH();
62        Int_t ntracks = TH->GetEntries();
63 //
64 //   Loop over tracks
65 //
66
67        Float_t pups[4]={0,0,0,0};
68        for (Int_t track=0; track<ntracks;track++) {
69            Int_t Nelt=0;           
70            gAlice->ResetHits();
71            Int_t nbytes += TH->GetEvent(track);
72
73            
74            if (MUON)  {
75                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
76                    mHit;
77                    mHit=(AliMUONHit*)MUON->NextHit()) 
78                {
79                    Int_t   nch   = mHit->Chamber();  // chamber number
80                    Float_t x     = mHit->X();        // x-pos of hit
81                    Float_t y     = mHit->Y();        // y-pos
82                    Float_t z     = mHit->Z();        // y-pos
83                    Float_t p  = mHit->Momentum();
84                    Float_t px = mHit->Px();
85                    Float_t py = mHit->Py();
86                    Float_t pz = mHit->Pz();
87                    
88                    if (nch != 1) continue;
89                    
90                    Int_t ipart = mHit->Particle();
91                    TParticle *Part;
92                    Int_t ftrack = mHit->Track();
93                    Part = gAlice->Particle(ftrack);
94                    Int_t ipart = Part->GetPdgCode();
95                    TParticle *Mother;
96                    Float_t px0 = Part->Px();
97                    Float_t py0 = Part->Py();               
98                    Float_t pz0 = Part->Pz();
99                    Float_t p0  = Part->P();
100                    
101                    Float_t e0=Part->Energy();
102                    
103                    if (ipart == kMuonPlus || ipart == kMuonMinus) {
104 //
105 // Branson Correction
106 //
107                        Float_t zch = z;
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 pzB = p * zb/TMath::Sqrt(zb*zb+yb*yb+xb*xb);
120                        Float_t pxB = pzB * xb/zb;
121                        Float_t pyB = pzB * yb/zb;
122                        Float_t theta 
123                            = TMath::ATan2(TMath::Sqrt(pxB*pxB+pyB*pyB), pzB) *
124                            180./TMath::Pi();
125                        
126                        
127 //
128 // Energy Correction
129 //
130 //
131                        
132                        Float_t corr=(p+CorrectP(p, theta))/p;
133                        pups[0]+=p*corr;
134                        pups[1]+=pxB*corr;
135                        pups[2]+=pyB*corr;
136                        pups[3]+=pzB*corr;
137                    }
138                } // hits 
139            } // if MUON 
140        } // tracks
141        Float_t mass =
142            TMath::Sqrt(pups[0]*pups[0]
143                        -pups[1]*pups[1]-pups[2]*pups[2]-pups[3]*pups[3]);
144        mass1->Fill(mass, 1.);
145        mass2->Fill(mass, 1.);
146        mass3->Fill(mass, 1.);          
147     } // event
148 //Create a canvas, set the view range, show histograms
149    TCanvas *c1 = new TCanvas("c1","Mass Spectra",400,10,600,700);
150
151    mass2->SetFillColor(42);
152    mass2->SetXTitle("Mass (GeV)");
153    mass2->Draw();
154 }
155
156
157
158 Float_t CorrectP(Float_t p, Float_t theta)
159 {
160     Float_t c;
161     
162     if (theta<3.) {
163 //W
164         if (p<15) {
165             c = 2.916+0.0156*p-0.00000866*p*p;
166         } else {
167             c = 2.998+0.0137*p;
168         }
169     } else {
170 //Cu+Fe+Cc+C
171         if (p<15) {
172             c = 2.464+0.00877*p-0.0000106*p*p;
173         } else {
174             c = 2.496+0.00817*p;
175         }
176     }
177 //
178 //  
179     c*=0.95;
180 //
181     return c;
182 }
183
184
185
186
187
188