Formatting changes.
[u/mrichter/AliRoot.git] / MUON / MUONacc.C
1 Bool_t SelectThis(AliMUONHit* mHit);
2
3 Bool_t SelectThis(Float_t x, Float_t y, Float_t d)
4 {
5 //    Float_t x     =    mHit->X();        // x-pos 
6 //    Float_t y     =    mHit->Y();        // y-pos
7     if (TMath::Abs(x) > d/2. && TMath::Abs(y) > d/2.) {
8         return 1;
9     } else {
10         return 0;
11     }
12 }
13
14 Bool_t SelectThat(Float_t x, Float_t y, Float_t d)
15 {
16     const Float_t phi0=60.*TMath::Pi()/180.;
17     
18     if (y<0.) return 1;
19     if (x<0.) x=-x;
20     
21     Float_t phi = TMath::ATan2(y,x);
22     Float_t dist = TMath::Sin(phi0-phi)*TMath::Sqrt(x*x+y*y);
23     if (TMath::Abs(dist) > d) {
24         return 1;
25     } else {
26         return 0 ;
27     }
28 }
29
30 void MUONacc (Float_t d=3., Int_t evNumber1=0, Int_t evNumber2=0) 
31 {
32 // Dynamically link some shared libs                    
33     Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
34                         56.1, 56.1, 59.6, 59.6};
35     Float_t rmax[14] = {91.5, 91.5, 122.5, 122.5, 158.3, 158.3, 260.0, 
36                         260.0, 260.0, 260.0, 293., 293, 293., 293.};
37     
38
39     if (gClassTable->GetID("AliRun") < 0) {
40         gROOT->LoadMacro("loadlibs.C");
41         loadlibs();
42     }
43 // Connect the Root Galice file containing Geometry, Kine and Hits
44
45     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
46
47     
48     if (!file) {
49         printf("\n Creating galice.root \n");
50         file = new TFile("galice.root");
51     } else {
52         printf("\n galice.root found in file list");
53     }
54     file->ls();
55 //
56     TDatabasePDG*  DataBase = new TDatabasePDG();
57     
58 // Get AliRun object from file or create it if not on file
59     if (!gAlice) {
60         gAlice = (AliRun*)(file->Get("gAlice"));
61         if (gAlice) printf("AliRun object found on file\n");
62         if (!gAlice) {
63             printf("\n create new gAlice object");
64             gAlice = new AliRun("gAlice","Alice test program");
65         }
66     }
67     
68 //  Create some histograms
69
70     TH1F *theta   =  new TH1F("theta","Theta distribution",180,0,180);
71     TH1F *emult   =  new TH1F("emult","Event Multiplicity",100,0,1000);   
72     TH2F *supp    =  new TH2F("supp", "Diselected",100, -200., 200, 100, -200., 200.);   
73
74     AliMUONChamber*  iChamber;
75     AliSegmentation* seg;
76     
77 //
78 //   Loop over events 
79 //
80     
81     for (Int_t nev=0; nev<= evNumber2; nev++) {
82         Int_t nparticles = gAlice->GetEvent(nev);
83         if (nev < evNumber1) continue;
84         if (nparticles <= 0) return;
85         
86         AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
87         
88         
89         TTree *TH = gAlice->TreeH();
90         Int_t ntracks = TH->GetEntries();
91         TClonesArray *fPartArray = gAlice->Particles();       
92         Int_t npart = fPartArray->GetEntriesFast();
93         
94         Int_t EvMult=0;
95         Int_t Nsel=0;
96 //
97 // Loop over primary particles (jpsi. upsilon, ...)
98 //       
99         Int_t prim=0;
100         
101         for (Int_t part=0; part<30000; part+=3) {
102             TParticle *MPart = gAlice->Particle(part);
103             Int_t mpart  = MPart->GetPdgCode();
104             Int_t child1 = MPart->GetFirstDaughter();
105             Int_t child2 = MPart->GetLastDaughter();    
106             Int_t mother = MPart->GetFirstMother();         
107 //
108 // Loop over muons
109 //
110             Int_t muon1=2*prim;
111             Int_t muon2=muon1+1;
112
113 //          printf("\n PDG Code %d %d %d %d %d %d %d %d \n", 
114 //                 part, mpart, child1, child2, mother,prim, muon1,muon2);
115
116
117             Bool_t selected[2]={1,1};
118             Int_t muons=0;
119             for (Int_t track=muon1; track<=muon2;track++) {
120                 gAlice->ResetHits();
121                 Int_t nbytes += TH->GetEvent(track);
122                 
123                 if (MUON)  {
124 //
125 //  Loop over hits
126 //
127                     for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
128                         mHit;
129                         mHit=(AliMUONHit*)MUON->NextHit()) 
130                     {
131                         Int_t   nch   =    mHit->Chamber();  // chamber number
132                         Float_t x     =    mHit->X();        // x-pos of hit
133                         Float_t y     =    mHit->Y();        // y-pos
134                         Float_t z     =    mHit->Z();        // z-pos
135                         Float_t Eloss =    mHit->Eloss();    // energy loss
136                         Float_t Theta =    mHit->Theta();    // theta
137                         Float_t Particle = mHit->Particle(); // Particle type
138                         Int_t itrack   = Int_t(mHit->Track());
139
140                         TParticle *thePart = gAlice->Particle(itrack);
141                         Float_t pTheta=thePart->Theta();
142                         
143                         if (Particle != kMuonPlus && Particle != kMuonMinus) continue;
144                         
145                         Float_t P   = mhit->Momentum();
146                         Float_t R   = TMath::Sqrt(x*x+y*y);
147                         TParticlePDG* Part = DataBase->GetParticle(Particle);
148                         Double_t mass = Part->Mass();
149                         theta->Fill(pTheta*180./TMath::Pi(),1.);
150                         
151                         if (nch>14) continue;
152 //                      printf("\n %f %f %f %d\n ", R, rmin[nch-1], rmax[nch-1], nch);
153 // Geometrical acceptance of frames                     
154 //                      selected[muons] = selected[muons]&&SelectThis(x,y,d);
155 // Cut on default rmin and rmax
156 //                      selected[muons] = selected[muons] &&( R > rmin[nch-1] && R < rmax[nch-1]);
157 // Cut on rmin
158 //                      selected[muons] = selected[muons] &&( R > rmin[nch-1]);
159 // Cut on rmax
160 //                      selected[muons] = selected[muons] &&( R < rmax[nch-1]);
161 //                      selected[muons] = selected[muons] &&( R > z*TMath::Tan(2.0*TMath::Pi()/180.));
162 //                      Bool_t ok = (z<970 && ( R > z*TMath::Tan(2.0*TMath::Pi()/180.))) ||
163 //                          (z>970 && ( R > 970 * TMath::Tan(2.0*TMath::Pi()/180.)
164 //                                      +(z-970)* TMath::Tan(1.6*TMath::Pi()/180.)));
165 //
166 //                      Bool_t ok = (z<970 && ( R < z*TMath::Tan(9.0*TMath::Pi()/180.))) ||
167 //                          (z>970 && ( R < 970 * TMath::Tan(9.0*TMath::Pi()/180.)
168 //                                      +(z-970)* TMath::Tan(12.0*TMath::Pi()/180.)));
169                         if (nch == 4 ||  nch == 7) {
170                             Bool_t ok = SelectThat(x,y,d);
171                         } else {
172                             Bool_t ok = kTRUE;
173                         }
174                         
175                         
176                             selected[muons] = selected[muons] && ok;
177                         if (!ok) {
178                             supp->Fill(x,y,1.);
179                             
180                         }
181                         
182
183                     } // hit loop
184                 } // if MUON
185                 muons++;
186             } // track loop
187             if (selected[0] && selected[1]) Nsel++;
188             prim++;
189         } // primary loop
190         printf("\n Selected %d %d \n", Nsel, prim);
191     }
192     
193 //Create a canvas, set the view range, show histograms
194     Int_t k;
195     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
196     c1->Divide(2,2);
197     c1->cd(1);
198     theta->Draw();
199     c1->cd(2);
200     supp->Draw();
201     
202 }
203
204
205
206
207
208
209
210
211