fe4da5cc |
1 | //////////////////////////////////////////////// |
2 | // Manager and hits classes for set:MUON // |
3 | //////////////////////////////////////////////// |
4 | |
5 | #include <TTUBE.h> |
a897a37a |
6 | #include <TBRIK.h> |
7 | #include <TRotMatrix.h> |
fe4da5cc |
8 | #include <TNode.h> |
a897a37a |
9 | #include <TTree.h> |
fe4da5cc |
10 | #include <TRandom.h> |
11 | #include <TObject.h> |
12 | #include <TVector.h> |
13 | #include <TObjArray.h> |
a897a37a |
14 | #include <TMinuit.h> |
15 | #include <TParticle.h> |
16 | #include <TROOT.h> |
17 | #include <TFile.h> |
18 | #include <TNtuple.h> |
19 | #include <TCanvas.h> |
20 | #include <TPad.h> |
21 | #include <TDirectory.h> |
22 | #include <TObjectTable.h> |
23 | #include <AliPDG.h> |
fe4da5cc |
24 | |
25 | #include "AliMUON.h" |
a897a37a |
26 | #include "TTUBE.h" |
27 | #include "AliMUONClusterFinder.h" |
fe4da5cc |
28 | #include "AliRun.h" |
29 | #include "AliMC.h" |
30 | #include "iostream.h" |
31 | #include "AliCallf77.h" |
32 | |
a897a37a |
33 | #ifndef WIN32 |
34 | # define reco_init reco_init_ |
35 | # define cutpxz cutpxz_ |
36 | # define sigmacut sigmacut_ |
37 | # define xpreci xpreci_ |
38 | # define ypreci ypreci_ |
39 | # define reconstmuon reconstmuon_ |
40 | # define trackf_read_geant trackf_read_geant_ |
41 | # define trackf_read_spoint trackf_read_spoint_ |
42 | # define chfill chfill_ |
43 | # define chfill2 chfill2_ |
44 | # define chf1 chf1_ |
45 | # define chfnt chfnt_ |
46 | # define hist_create hist_create_ |
47 | # define hist_closed hist_closed_ |
48 | # define rndm rndm_ |
49 | # define fcn fcn_ |
50 | # define trackf_fit trackf_fit_ |
51 | # define prec_fit prec_fit_ |
52 | # define fcnfit fcnfit_ |
53 | # define reco_term reco_term_ |
54 | #else |
55 | # define reco_init RECO_INIT |
56 | # define cutpxz CUTPXZ |
57 | # define sigmacut SIGMACUT |
58 | # define xpreci XPRECI |
59 | # define ypreci YPRECI |
60 | # define reconstmuon RECONSTMUON |
61 | # define trackf_read_geant TRACKF_READ_GEANT |
62 | # define trackf_read_spoint TRACKF_READ_SPOINT |
63 | # define chfill CHFILL |
64 | # define chfill2 CHFILL2 |
65 | # define chf1 CHF1 |
66 | # define chfnt CHFNT |
67 | # define hist_create HIST_CREATE |
68 | # define hist_closed HIST_CLOSED |
69 | # define rndm RNDM |
70 | # define fcn FCN |
71 | # define trackf_fit TRACKF_FIT |
72 | # define prec_fit PREC_FIT |
73 | # define fcnfit FCNFIT |
74 | # define reco_term RECO_TERM |
75 | #endif |
76 | |
77 | extern "C" |
78 | { |
79 | void type_of_call reco_init(Double_t &, Double_t &, Double_t &); |
80 | void type_of_call reco_term(); |
81 | void type_of_call cutpxz(Double_t &); |
82 | void type_of_call sigmacut(Double_t &); |
83 | void type_of_call xpreci(Double_t &); |
84 | void type_of_call ypreci(Double_t &); |
85 | void type_of_call reconstmuon(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &); |
86 | void type_of_call trackf_read_geant(Int_t *, Double_t *, Double_t *, Double_t *, Int_t *, Int_t *, Double_t *, Double_t *, Double_t *, Double_t *,Int_t &, Double_t *, Double_t *, Double_t *, Int_t &, Int_t &, Double_t *, Double_t *, Double_t *, Double_t *); |
87 | void type_of_call trackf_read_spoint(Int_t *, Double_t *, Double_t *, Double_t *, Int_t *, Int_t *, Double_t *, Double_t *, Double_t *, Double_t *,Int_t &, Double_t *, Double_t *, Double_t *, Int_t &, Int_t &, Double_t *, Double_t *, Double_t *, Double_t *); |
88 | void type_of_call chfill(Int_t &, Float_t &, Float_t &, Float_t &); |
89 | void type_of_call chfill2(Int_t &, Float_t &, Float_t &, Float_t &); |
90 | void type_of_call chf1(Int_t &, Float_t &, Float_t &); |
91 | void type_of_call chfnt(Int_t &, Int_t &, Int_t *, Int_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *); |
92 | void type_of_call hist_create(); |
93 | void type_of_call hist_closed(); |
94 | void type_of_call fcnf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t); |
95 | void type_of_call fcn(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &); |
96 | void type_of_call trackf_fit(Int_t &, Double_t *, Double_t *, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &); |
97 | void type_of_call prec_fit(Double_t &, Double_t &, Double_t &, Double_t &, Double_t&, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &); |
98 | void type_of_call fcnfitf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t); |
99 | void type_of_call fcnfit(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &); |
100 | Float_t type_of_call rndm() {return gRandom->Rndm();} |
101 | } |
102 | |
fe4da5cc |
103 | // Static variables for the pad-hit iterator routines |
104 | static Int_t sMaxIterPad=0; |
105 | static Int_t sCurIterPad=0; |
a6f39961 |
106 | static TTree *TrH1; |
a897a37a |
107 | static TTree *TK1; |
108 | static TClonesArray *fHits2; //Listof hits for one track only |
109 | static TClonesArray *fClusters2; //List of clusters for one track only |
110 | static TClonesArray *fParticles2; //List of particles in the Kine tree |
fe4da5cc |
111 | ClassImp(AliMUON) |
fe4da5cc |
112 | //___________________________________________ |
113 | AliMUON::AliMUON() |
114 | { |
115 | fIshunt = 0; |
116 | fHits = 0; |
117 | fClusters = 0; |
118 | fNclusters = 0; |
119 | fDchambers = 0; |
fe4da5cc |
120 | fNdch = 0; |
a897a37a |
121 | fRawClusters= 0; |
122 | fNrawch = 0; |
123 | fCathCorrel= 0; |
124 | fNcorch = 0; |
125 | fTreeC = 0; |
126 | |
127 | // modifs perso |
128 | fSPxzCut = 0; |
129 | fSSigmaCut = 0; |
130 | fSXPrec = 0; |
131 | fSYPrec = 0; |
fe4da5cc |
132 | } |
133 | |
134 | //___________________________________________ |
135 | AliMUON::AliMUON(const char *name, const char *title) |
136 | : AliDetector(name,title) |
137 | { |
138 | //Begin_Html |
139 | /* |
a897a37a |
140 | <img src="gif/alimuon.gif"> |
fe4da5cc |
141 | */ |
142 | //End_Html |
143 | |
a897a37a |
144 | fHits = new TClonesArray("AliMUONhit",1000); |
fe4da5cc |
145 | fClusters = new TClonesArray("AliMUONcluster",10000); |
146 | fNclusters = 0; |
147 | fIshunt = 0; |
148 | |
a897a37a |
149 | fNdch = new Int_t[10]; |
fe4da5cc |
150 | |
a897a37a |
151 | fDchambers = new TObjArray(10); |
fe4da5cc |
152 | |
153 | Int_t i; |
154 | |
a897a37a |
155 | for (i=0; i<10 ;i++) { |
fe4da5cc |
156 | (*fDchambers)[i] = new TClonesArray("AliMUONdigit",10000); |
157 | fNdch[i]=0; |
158 | } |
159 | |
a897a37a |
160 | fNrawch = new Int_t[10]; |
161 | |
162 | fRawClusters = new TObjArray(10); |
163 | |
164 | for (i=0; i<10 ;i++) { |
165 | (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000); |
166 | fNrawch[i]=0; |
167 | } |
168 | |
169 | fNcorch = new Int_t[10]; |
170 | fCathCorrel = new TObjArray(10); |
171 | for (i=0; i<10 ;i++) { |
172 | (*fCathCorrel)[i] = new TClonesArray("AliMUONcorrelation",1000); |
173 | fNcorch[i]=0; |
174 | } |
175 | |
176 | fTreeC = 0; |
fe4da5cc |
177 | |
178 | // |
179 | // Transport angular cut |
180 | fAccCut=0; |
181 | fAccMin=2; |
182 | fAccMax=9; |
183 | |
a897a37a |
184 | // modifs perso |
185 | fSPxzCut = 3.0; |
186 | fSSigmaCut = 1.0; |
187 | fSXPrec = 0.01; |
188 | fSYPrec = 0.144; |
e3a4d40e |
189 | |
fe4da5cc |
190 | SetMarkerColor(kRed); |
191 | } |
192 | |
193 | //___________________________________________ |
194 | AliMUON::~AliMUON() |
195 | { |
a897a37a |
196 | |
197 | printf("Calling AliMUON destructor !!!\n"); |
198 | |
a6f39961 |
199 | Int_t i; |
fe4da5cc |
200 | fIshunt = 0; |
201 | delete fHits; |
202 | delete fClusters; |
a897a37a |
203 | delete fTreeC; |
fe4da5cc |
204 | |
a6f39961 |
205 | for (i=0;i<10;i++) { |
a897a37a |
206 | delete (*fDchambers)[i]; |
207 | fNdch[i]=0; |
208 | } |
209 | delete fDchambers; |
210 | |
a6f39961 |
211 | for (i=0;i<10;i++) { |
a897a37a |
212 | delete (*fRawClusters)[i]; |
213 | fNrawch[i]=0; |
214 | } |
215 | delete fRawClusters; |
216 | |
a6f39961 |
217 | for (i=0;i<10;i++) { |
a897a37a |
218 | delete (*fCathCorrel)[i]; |
219 | fNcorch[i]=0; |
220 | } |
221 | delete fCathCorrel; |
fe4da5cc |
222 | } |
223 | |
224 | //___________________________________________ |
225 | void AliMUON::AddHit(Int_t track, Int_t *vol, Float_t *hits) |
226 | { |
227 | TClonesArray &lhits = *fHits; |
228 | new(lhits[fNhits++]) AliMUONhit(fIshunt,track,vol,hits); |
229 | } |
230 | //___________________________________________ |
231 | void AliMUON::AddCluster(Int_t *clhits) |
232 | { |
233 | TClonesArray &lclusters = *fClusters; |
234 | new(lclusters[fNclusters++]) AliMUONcluster(clhits); |
235 | } |
236 | //_____________________________________________________________________________ |
237 | void AliMUON::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits) |
238 | { |
239 | // |
240 | // Add a MUON digit to the list |
241 | // |
242 | |
243 | TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]); |
244 | new(ldigits[fNdch[id]++]) AliMUONdigit(tracks,charges,digits); |
245 | } |
246 | |
a897a37a |
247 | //_____________________________________________________________________________ |
248 | void AliMUON::AddRawCluster(Int_t id, const AliMUONRawCluster& c) |
249 | { |
250 | // |
251 | // Add a MUON digit to the list |
252 | // |
253 | |
254 | TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]); |
255 | new(lrawcl[fNrawch[id]++]) AliMUONRawCluster(c); |
256 | } |
257 | //_____________________________________________________________________________ |
258 | void AliMUON::AddCathCorrel(Int_t id, Int_t *idx, Float_t *x, Float_t *y) |
259 | { |
260 | // |
261 | // Add a MUON digit to the list |
262 | // |
263 | |
264 | TClonesArray &lcorrel = *((TClonesArray*)(*fCathCorrel)[id]); |
265 | new(lcorrel[fNcorch[id]++]) AliMUONcorrelation(idx,x,y); |
266 | } |
267 | |
fe4da5cc |
268 | //___________________________________________ |
269 | void AliMUON::BuildGeometry() |
270 | { |
a897a37a |
271 | TNode *Node, *NodeF, *Top; |
272 | const int kColorMUON = kBlue; |
273 | // |
274 | Top=gAlice->GetGeometry()->GetNode("alice"); |
275 | // MUON |
276 | // |
277 | // z-Positions of Chambers |
278 | const Float_t cz[5]={511., 686., 971., 1245., 1445.}; |
279 | // |
280 | // inner diameter |
281 | const Float_t dmi[5]={ 35., 47., 67., 86., 100.}; |
282 | // |
283 | // outer diameter |
284 | const Float_t dma[5]={183., 245., 346., 520., 520.}; |
285 | |
286 | TRotMatrix* rot000 = new TRotMatrix("Rot000"," ", 90, 0, 90, 90, 0, 0); |
287 | TRotMatrix* rot090 = new TRotMatrix("Rot090"," ", 90, 90, 90,180, 0, 0); |
288 | TRotMatrix* rot180 = new TRotMatrix("Rot180"," ", 90,180, 90,270, 0, 0); |
289 | TRotMatrix* rot270 = new TRotMatrix("Rot270"," ", 90,270, 90, 0, 0, 0); |
290 | |
291 | |
292 | float rmin, rmax, dx, dy, dz, dr, zpos; |
293 | float dzc=4.; |
294 | char NameChamber[9], NameSense[9], NameFrame[9], NameNode[7]; |
295 | for (Int_t i=0; i<5; i++) { |
296 | for (Int_t j=0; j<2; j++) { |
297 | Int_t id=2*i+j+1; |
298 | if (j==0) { |
299 | zpos=cz[i]-dzc; |
300 | } else { |
301 | zpos=cz[i]+dzc; |
302 | } |
303 | |
304 | |
305 | sprintf(NameChamber,"C_MUON%d",id); |
306 | sprintf(NameSense,"S_MUON%d",id); |
307 | sprintf(NameFrame,"F_MUON%d",id); |
308 | rmin = dmi[i]/2.-3; |
309 | rmax = dma[i]/2.+3; |
310 | new TTUBE(NameChamber,"Mother","void",rmin,rmax,0.25,1.); |
311 | rmin = dmi[i]/2.; |
312 | rmax = dma[i]/2.; |
313 | new TTUBE(NameSense,"Sens. region","void",rmin,rmax,0.25, 1.); |
314 | dx=(rmax-rmin)/2; |
315 | dy=3.; |
316 | dz=0.25; |
317 | TBRIK* FMUON = new TBRIK(NameFrame,"Frame","void",dx,dy,dz); |
318 | Top->cd(); |
319 | sprintf(NameNode,"MUON%d",100+id); |
320 | Node = new TNode(NameNode,"ChamberNode",NameChamber,0,0,zpos,""); |
321 | Node->SetLineColor(kColorMUON); |
322 | fNodes->Add(Node); |
323 | Node->cd(); |
324 | sprintf(NameNode,"MUON%d",200+id); |
325 | Node = new TNode(NameNode,"Sens. Region Node",NameSense,0,0,0,""); |
326 | Node->SetLineColor(kColorMUON); |
b0236364 |
327 | fNodes->Add(Node); |
a897a37a |
328 | Node->cd(); |
329 | dr=dx+rmin; |
330 | sprintf(NameNode,"MUON%d",300+id); |
331 | NodeF = new TNode(NameNode,"Frame0",FMUON,dr, 0, 0,rot000,""); |
332 | NodeF->SetLineColor(kColorMUON); |
b0236364 |
333 | fNodes->Add(NodeF); |
a897a37a |
334 | Node->cd(); |
335 | sprintf(NameNode,"MUON%d",400+id); |
336 | NodeF = new TNode(NameNode,"Frame1",FMUON,0 ,dr,0,rot090,""); |
337 | NodeF->SetLineColor(kColorMUON); |
b0236364 |
338 | fNodes->Add(NodeF); |
a897a37a |
339 | Node->cd(); |
340 | sprintf(NameNode,"MUON%d",500+id); |
341 | NodeF = new TNode(NameNode,"Frame2",FMUON,-dr,0,0,rot180,""); |
342 | NodeF->SetLineColor(kColorMUON); |
b0236364 |
343 | fNodes->Add(NodeF); |
a897a37a |
344 | Node ->cd(); |
345 | sprintf(NameNode,"MUON%d",600+id); |
346 | NodeF = new TNode(NameNode,"Frame3",FMUON,0,-dr,0,rot270,""); |
347 | NodeF->SetLineColor(kColorMUON); |
b0236364 |
348 | fNodes->Add(NodeF); |
a897a37a |
349 | } |
350 | } |
fe4da5cc |
351 | } |
352 | |
a897a37a |
353 | |
fe4da5cc |
354 | //___________________________________________ |
355 | Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t ) |
356 | { |
357 | return 9999; |
358 | } |
359 | |
360 | //___________________________________________ |
361 | void AliMUON::MakeBranch(Option_t* option) |
362 | { |
363 | // Create Tree branches for the MUON. |
364 | |
365 | const Int_t buffersize = 4000; |
a897a37a |
366 | char branchname[30]; |
fe4da5cc |
367 | sprintf(branchname,"%sCluster",GetName()); |
368 | |
369 | AliDetector::MakeBranch(option); |
370 | |
371 | if (fClusters && gAlice->TreeH()) { |
372 | gAlice->TreeH()->Branch(branchname,&fClusters, buffersize); |
373 | printf("Making Branch %s for clusters\n",branchname); |
374 | } |
375 | |
376 | // one branch for digits per chamber |
377 | Int_t i; |
378 | |
379 | for (i=0; i<10 ;i++) { |
380 | sprintf(branchname,"%sDigits%d",GetName(),i+1); |
381 | |
382 | if (fDchambers && gAlice->TreeD()) { |
383 | gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), buffersize); |
384 | printf("Making Branch %s for digits in chamber %d\n",branchname,i+1); |
385 | } |
386 | } |
a897a37a |
387 | |
e3a4d40e |
388 | //printf("Make Branch - TreeR address %p\n",gAlice->TreeR()); |
a897a37a |
389 | |
390 | // one branch for raw clusters per chamber |
391 | for (i=0; i<10 ;i++) { |
392 | sprintf(branchname,"%sRawClusters%d",GetName(),i+1); |
393 | |
394 | if (fRawClusters && gAlice->TreeR()) { |
395 | gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), buffersize); |
396 | printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1); |
397 | } |
398 | } |
399 | |
fe4da5cc |
400 | } |
401 | |
402 | //___________________________________________ |
403 | void AliMUON::SetTreeAddress() |
404 | { |
405 | // Set branch address for the Hits and Digits Tree. |
a897a37a |
406 | char branchname[30]; |
fe4da5cc |
407 | AliDetector::SetTreeAddress(); |
408 | |
409 | TBranch *branch; |
410 | TTree *treeH = gAlice->TreeH(); |
411 | TTree *treeD = gAlice->TreeD(); |
a897a37a |
412 | TTree *treeR = gAlice->TreeR(); |
fe4da5cc |
413 | |
414 | if (treeH) { |
415 | if (fClusters) { |
416 | branch = treeH->GetBranch("MUONCluster"); |
417 | if (branch) branch->SetAddress(&fClusters); |
418 | } |
419 | } |
420 | |
421 | if (treeD) { |
422 | for (int i=0; i<10; i++) { |
423 | sprintf(branchname,"%sDigits%d",GetName(),i+1); |
424 | if (fDchambers) { |
425 | branch = treeD->GetBranch(branchname); |
426 | if (branch) branch->SetAddress(&((*fDchambers)[i])); |
427 | } |
428 | } |
429 | } |
a897a37a |
430 | |
431 | // printf("SetTreeAddress --- treeR address %p \n",treeR); |
432 | |
433 | if (treeR) { |
434 | for (int i=0; i<10; i++) { |
435 | sprintf(branchname,"%sRawClusters%d",GetName(),i+1); |
436 | if (fRawClusters) { |
437 | branch = treeR->GetBranch(branchname); |
438 | if (branch) branch->SetAddress(&((*fRawClusters)[i])); |
439 | } |
440 | } |
441 | } |
442 | |
fe4da5cc |
443 | } |
444 | //___________________________________________ |
445 | void AliMUON::ResetHits() |
446 | { |
447 | // Reset number of clusters and the cluster array for this detector |
448 | AliDetector::ResetHits(); |
449 | fNclusters = 0; |
450 | if (fClusters) fClusters->Clear(); |
451 | } |
452 | |
453 | //____________________________________________ |
454 | void AliMUON::ResetDigits() |
455 | { |
456 | // |
457 | // Reset number of digits and the digits array for this detector |
458 | // |
459 | for ( int i=0;i<10;i++ ) { |
a897a37a |
460 | if ((*fDchambers)[i]) ((TClonesArray*)(*fDchambers)[i])->Clear(); |
fe4da5cc |
461 | if (fNdch) fNdch[i]=0; |
462 | } |
463 | } |
a897a37a |
464 | //____________________________________________ |
465 | void AliMUON::ResetRawClusters() |
466 | { |
467 | // |
468 | // Reset number of raw clusters and the raw clust array for this detector |
469 | // |
470 | for ( int i=0;i<10;i++ ) { |
471 | if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear(); |
472 | if (fNrawch) fNrawch[i]=0; |
473 | } |
474 | } |
475 | //____________________________________________ |
476 | void AliMUON::ResetCorrelation() |
477 | { |
478 | // |
479 | // Reset number of correl clusters and the correl clust array for |
480 | // this detector |
481 | // |
482 | for ( int i=0;i<10;i++ ) { |
483 | if ((*fCathCorrel)[i]) ((TClonesArray*)(*fCathCorrel)[i])->Clear(); |
484 | if (fNcorch) fNcorch[i]=0; |
485 | } |
486 | } |
487 | |
fe4da5cc |
488 | //___________________________________________ |
489 | |
490 | void AliMUON::SetPADSIZ(Int_t id, Int_t isec, Float_t p1, Float_t p2) |
491 | { |
492 | Int_t i=2*(id-1); |
493 | ((AliMUONchamber*) (*fChambers)[i]) ->SetPADSIZ(isec,p1,p2); |
494 | ((AliMUONchamber*) (*fChambers)[i+1])->SetPADSIZ(isec,p1,p2); |
495 | } |
496 | |
497 | //___________________________________________ |
a897a37a |
498 | void AliMUON::SetChargeSlope(Int_t id, Float_t p1) |
fe4da5cc |
499 | { |
500 | Int_t i=2*(id-1); |
a897a37a |
501 | ((AliMUONchamber*) (*fChambers)[i])->SetChargeSlope(p1); |
502 | ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSlope(p1); |
fe4da5cc |
503 | } |
504 | |
505 | //___________________________________________ |
a897a37a |
506 | void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2) |
fe4da5cc |
507 | { |
508 | Int_t i=2*(id-1); |
a897a37a |
509 | ((AliMUONchamber*) (*fChambers)[i])->SetChargeSpread(p1,p2); |
510 | ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSpread(p1,p2); |
fe4da5cc |
511 | } |
512 | |
513 | //___________________________________________ |
a897a37a |
514 | void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1) |
fe4da5cc |
515 | { |
516 | Int_t i=2*(id-1); |
a897a37a |
517 | ((AliMUONchamber*) (*fChambers)[i])->SetSigmaIntegration(p1); |
518 | ((AliMUONchamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1); |
fe4da5cc |
519 | } |
520 | |
521 | //___________________________________________ |
a897a37a |
522 | void AliMUON::SetMaxAdc(Int_t id, Float_t p1) |
fe4da5cc |
523 | { |
524 | Int_t i=2*(id-1); |
a897a37a |
525 | ((AliMUONchamber*) (*fChambers)[i])->SetMaxAdc(p1); |
526 | ((AliMUONchamber*) (*fChambers)[i+1])->SetMaxAdc(p1); |
fe4da5cc |
527 | } |
528 | |
529 | //___________________________________________ |
a897a37a |
530 | void AliMUON::SetMaxStepGas(Float_t p1) |
fe4da5cc |
531 | { |
532 | fMaxStepGas=p1; |
533 | } |
534 | |
535 | //___________________________________________ |
a897a37a |
536 | void AliMUON::SetMaxStepAlu(Float_t p1) |
fe4da5cc |
537 | { |
538 | fMaxStepAlu=p1; |
539 | } |
540 | |
541 | //___________________________________________ |
a897a37a |
542 | void AliMUON::SetMaxDestepGas(Float_t p1) |
fe4da5cc |
543 | { |
544 | fMaxDestepGas=p1; |
545 | } |
546 | |
547 | //___________________________________________ |
a897a37a |
548 | void AliMUON::SetMaxDestepAlu(Float_t p1) |
fe4da5cc |
549 | { |
550 | fMaxDestepAlu=p1; |
551 | } |
552 | //___________________________________________ |
a897a37a |
553 | void AliMUON::SetMuonAcc(Bool_t acc, Float_t angmin, Float_t angmax) |
fe4da5cc |
554 | { |
555 | fAccCut=acc; |
556 | fAccMin=angmin; |
557 | fAccMax=angmax; |
558 | } |
559 | //___________________________________________ |
560 | void AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliMUONsegmentation *segmentation) |
561 | { |
562 | ((AliMUONchamber*) (*fChambers)[id])->SegmentationModel(isec, segmentation); |
563 | |
564 | } |
565 | //___________________________________________ |
566 | void AliMUON::SetResponseModel(Int_t id, AliMUONresponse *response) |
567 | { |
568 | ((AliMUONchamber*) (*fChambers)[id])->ResponseModel(response); |
569 | } |
570 | |
a897a37a |
571 | void AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinder *reconst) |
572 | { |
573 | ((AliMUONchamber*) (*fChambers)[id])->ReconstructionModel(reconst); |
574 | } |
575 | |
fe4da5cc |
576 | void AliMUON::SetNsec(Int_t id, Int_t nsec) |
577 | { |
578 | ((AliMUONchamber*) (*fChambers)[id])->SetNsec(nsec); |
579 | } |
580 | |
581 | |
582 | //___________________________________________ |
583 | |
584 | void AliMUON::StepManager() |
585 | { |
586 | printf("Dummy version of muon step -- it should never happen!!\n"); |
e3a4d40e |
587 | /* |
fe4da5cc |
588 | const Float_t kRaddeg = 180/TMath::Pi(); |
a897a37a |
589 | AliMC* pMC = AliMC::GetMC(); |
fe4da5cc |
590 | Int_t nsec, ipart; |
591 | Float_t x[4], p[4]; |
a897a37a |
592 | Float_t pt, th0, th2; |
fe4da5cc |
593 | char proc[5]; |
594 | if(fAccCut) { |
a897a37a |
595 | if((nsec=pMC->NSecondaries())>0) { |
596 | pMC->ProdProcess(proc); |
597 | if((pMC->TrackPid()==443 || pMC->TrackPid()==553) && !strcmp(proc,"DCAY")) { |
fe4da5cc |
598 | // |
599 | // Check angular acceptance |
e3a4d40e |
600 | // --- and have muons from resonance decays in the wanted window --- |
fe4da5cc |
601 | if(nsec != 2) { |
602 | printf(" AliMUON::StepManager: Strange resonance Decay into %d particles\n",nsec); |
a897a37a |
603 | pMC->StopEvent(); |
fe4da5cc |
604 | } else { |
a897a37a |
605 | pMC->GetSecondary(0,ipart,x,p); |
fe4da5cc |
606 | pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]); |
607 | th0 = TMath::ATan2(pt,p[2])*kRaddeg; |
e3a4d40e |
608 | pMC->GetSecondary(1,ipart,x,p); |
fe4da5cc |
609 | pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]); |
a897a37a |
610 | th2 = TMath::ATan2(pt,p[2])*kRaddeg; |
fe4da5cc |
611 | if(!(fAccMin < th0 && th0 < fAccMax) || |
a897a37a |
612 | !(fAccMin < th2 && th2 < fAccMax)) |
613 | pMC->StopEvent(); |
fe4da5cc |
614 | } |
615 | } |
616 | } |
617 | } |
e3a4d40e |
618 | */ |
fe4da5cc |
619 | } |
a897a37a |
620 | |
621 | void AliMUON::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol) |
fe4da5cc |
622 | { |
623 | // |
a897a37a |
624 | // Calls the charge disintegration method of the current chamber and adds |
625 | // the simulated cluster to the root treee |
fe4da5cc |
626 | // |
a897a37a |
627 | Int_t clhits[7]; |
628 | Float_t newclust[6][500]; |
629 | Int_t nnew; |
630 | |
631 | |
fe4da5cc |
632 | // |
a897a37a |
633 | // Integrated pulse height on chamber |
634 | |
635 | |
636 | clhits[0]=fNhits+1; |
fe4da5cc |
637 | // |
a897a37a |
638 | // |
639 | ((AliMUONchamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust); |
640 | // printf("\n Add new clusters %d %f \n", nnew, eloss*1.e9); |
641 | Int_t ic=0; |
642 | |
643 | // |
644 | // Add new clusters |
645 | for (Int_t i=0; i<nnew; i++) { |
646 | if (Int_t(newclust[3][i]) > 0) { |
647 | ic++; |
648 | // Cathode plane |
649 | clhits[1] = Int_t(newclust[5][i]); |
650 | // Cluster Charge |
651 | clhits[2] = Int_t(newclust[0][i]); |
652 | // Pad: ix |
653 | clhits[3] = Int_t(newclust[1][i]); |
654 | // Pad: iy |
655 | clhits[4] = Int_t(newclust[2][i]); |
656 | // Pad: charge |
657 | clhits[5] = Int_t(newclust[3][i]); |
658 | // Pad: chamber sector |
659 | clhits[6] = Int_t(newclust[4][i]); |
660 | |
661 | AddCluster(clhits); |
662 | } |
663 | } |
664 | // printf("\n %d new clusters added", ic); |
665 | } |
666 | |
e3a4d40e |
667 | void AliMUON::Digitise(Int_t nev,Int_t bgr_ev,Option_t *option, Option_t *,Text_t *filename) |
a897a37a |
668 | { |
669 | // keep galice.root for signal and name differently the file for |
670 | // background when add! otherwise the track info for signal will be lost ! |
671 | |
a6f39961 |
672 | static Bool_t first=kTRUE; |
673 | // static TTree *TrH1; |
a897a37a |
674 | static TFile *File; |
675 | char *Add = strstr(option,"Add"); |
676 | //char *listoftracks = strstr(opt,"listoftracks"); |
677 | |
678 | AliMUONchamber* iChamber; |
679 | AliMUONsegmentation* segmentation; |
680 | |
681 | |
682 | Int_t trk[50]; |
683 | Int_t chtrk[50]; |
684 | TObjArray *list=new TObjArray; |
685 | static TClonesArray *p_adr=0; |
686 | if(!p_adr) p_adr=new TClonesArray("TVector",1000); |
687 | Int_t digits[5]; |
688 | |
689 | AliMUON *MUON = (AliMUON *) gAlice->GetModule("MUON"); |
690 | AliMUONHitMap * HitMap[10]; |
691 | for (Int_t i=0; i<10; i++) {HitMap[i]=0;} |
692 | if (Add ) { |
693 | if(first) { |
694 | fFileName=filename; |
695 | cout<<"filename"<<fFileName<<endl; |
696 | File=new TFile(fFileName); |
697 | cout<<"I have opened "<<fFileName<<" file "<<endl; |
698 | fHits2 = new TClonesArray("AliMUONhit",1000 ); |
699 | fClusters2 = new TClonesArray("AliMUONcluster",10000); |
700 | } |
a6f39961 |
701 | first=kFALSE; |
a897a37a |
702 | File->cd(); |
703 | //File->ls(); |
704 | // Get Hits Tree header from file |
705 | if(fHits2) fHits2->Clear(); |
706 | if(fClusters2) fClusters2->Clear(); |
a6f39961 |
707 | if(TrH1) delete TrH1; |
708 | TrH1=0; |
a897a37a |
709 | |
710 | char treeName[20]; |
711 | sprintf(treeName,"TreeH%d",bgr_ev); |
a6f39961 |
712 | TrH1 = (TTree*)gDirectory->Get(treeName); |
713 | //printf("TrH1 %p of treename %s for event %d \n",TrH1,treeName,bgr_ev); |
a897a37a |
714 | |
a6f39961 |
715 | if (!TrH1) { |
a897a37a |
716 | printf("ERROR: cannot find Hits Tree for event:%d\n",bgr_ev); |
717 | } |
718 | // Set branch addresses |
719 | TBranch *branch; |
720 | char branchname[20]; |
721 | sprintf(branchname,"%s",GetName()); |
a6f39961 |
722 | if (TrH1 && fHits2) { |
723 | branch = TrH1->GetBranch(branchname); |
a897a37a |
724 | if (branch) branch->SetAddress(&fHits2); |
725 | } |
a6f39961 |
726 | if (TrH1 && fClusters2) { |
727 | branch = TrH1->GetBranch("MUONCluster"); |
a897a37a |
728 | if (branch) branch->SetAddress(&fClusters2); |
729 | } |
730 | // test |
a6f39961 |
731 | //Int_t ntracks1 =(Int_t)TrH1->GetEntries(); |
a897a37a |
732 | //printf("background - ntracks1 - %d\n",ntracks1); |
733 | } |
734 | // |
735 | // loop over cathodes |
736 | // |
737 | AliMUONHitMap* hm; |
738 | Int_t countadr=0; |
739 | for (int icat=0; icat<2; icat++) { |
a897a37a |
740 | Int_t counter=0; |
741 | for (Int_t i =0; i<10; i++) { |
742 | iChamber=(AliMUONchamber*) (*fChambers)[i]; |
743 | if (iChamber->Nsec()==1 && icat==1) { |
fe4da5cc |
744 | continue; |
a897a37a |
745 | } else { |
746 | segmentation=iChamber->GetSegmentationModel(icat+1); |
fe4da5cc |
747 | } |
a897a37a |
748 | HitMap[i] = new AliMUONHitMapA1(segmentation, list); |
749 | } |
750 | //printf("Start loop over tracks \n"); |
751 | // |
752 | // Loop over tracks |
753 | // |
754 | |
755 | TTree *TH = gAlice->TreeH(); |
756 | Int_t ntracks =(Int_t) TH->GetEntries(); |
757 | //printf("signal - ntracks %d\n",ntracks); |
758 | Int_t nmuon[10]={0,0,0,0,0,0,0,0,0,0}; |
759 | Float_t xhit[10][2]; |
760 | Float_t yhit[10][2]; |
761 | |
762 | for (Int_t track=0; track<ntracks; track++) { |
763 | gAlice->ResetHits(); |
764 | TH->GetEvent(track); |
765 | |
766 | // |
767 | // Loop over hits |
768 | for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); |
769 | mHit; |
770 | mHit=(AliMUONhit*)MUON->NextHit()) |
771 | { |
772 | Int_t nch = mHit->fChamber-1; // chamber number |
773 | if (nch >9) continue; |
774 | iChamber = &(MUON->Chamber(nch)); |
775 | Int_t rmin = (Int_t)iChamber->RInner(); |
776 | Int_t rmax = (Int_t)iChamber->ROuter(); |
777 | // new 17.07.99 |
778 | if (Add) { |
779 | |
780 | if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) { |
781 | xhit[nch][nmuon[nch]]=mHit->fX; |
782 | yhit[nch][nmuon[nch]]=mHit->fY; |
783 | nmuon[nch]++; |
784 | if (nmuon[nch] >2) printf("nmuon %d\n",nmuon[nch]); |
785 | |
786 | } |
fe4da5cc |
787 | } |
a897a37a |
788 | |
789 | |
790 | |
791 | |
792 | // |
793 | // Loop over pad hits |
794 | for (AliMUONcluster* mPad= |
795 | (AliMUONcluster*)MUON->FirstPad(mHit,fClusters); |
796 | mPad; |
797 | mPad=(AliMUONcluster*)MUON->NextPad(fClusters)) |
798 | { |
799 | Int_t cathode = mPad->fCathode; // cathode number |
800 | Int_t ipx = mPad->fPadX; // pad number on X |
801 | Int_t ipy = mPad->fPadY; // pad number on Y |
802 | Int_t iqpad = Int_t(mPad->fQpad*kScale);// charge per pad |
803 | // Int_t iqpad = mPad->fQpad; // charge per pad |
804 | // |
805 | // |
806 | |
807 | if (cathode != (icat+1)) continue; |
808 | // fill the info array |
809 | Float_t thex, they; |
810 | segmentation=iChamber->GetSegmentationModel(cathode); |
811 | segmentation->GetPadCxy(ipx,ipy,thex,they); |
812 | Float_t rpad=TMath::Sqrt(thex*thex+they*they); |
813 | if (rpad < rmin || iqpad ==0 || rpad > rmax) continue; |
814 | |
815 | new((*p_adr)[countadr++]) TVector(2); |
816 | TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]); |
817 | trinfo(0)=(Float_t)track; |
818 | trinfo(1)=(Float_t)iqpad; |
819 | |
820 | digits[0]=ipx; |
821 | digits[1]=ipy; |
822 | digits[2]=iqpad; |
823 | digits[3]=iqpad; |
824 | if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) { |
825 | digits[4]=mPad->fHitNumber; |
826 | } else digits[4]=-1; |
827 | |
828 | AliMUONlist* pdigit; |
829 | // build the list of fired pads and update the info |
830 | if (!HitMap[nch]->TestHit(ipx, ipy)) { |
831 | |
832 | list->AddAtAndExpand( |
833 | new AliMUONlist(nch,digits),counter); |
834 | |
835 | HitMap[nch]->SetHit(ipx, ipy, counter); |
836 | counter++; |
837 | pdigit=(AliMUONlist*)list->At(list->GetLast()); |
838 | // list of tracks |
839 | TObjArray *trlist=(TObjArray*)pdigit->TrackList(); |
840 | trlist->Add(&trinfo); |
841 | } else { |
842 | pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy); |
843 | // update charge |
844 | (*pdigit).fSignal+=iqpad; |
845 | (*pdigit).fPhysics+=iqpad; |
846 | // update list of tracks |
847 | TObjArray* trlist=(TObjArray*)pdigit->TrackList(); |
848 | Int_t last_entry=trlist->GetLast(); |
849 | TVector *ptrk_p=(TVector*)trlist->At(last_entry); |
850 | TVector &ptrk=*ptrk_p; |
851 | Int_t last_track=Int_t(ptrk(0)); |
852 | Int_t last_charge=Int_t(ptrk(1)); |
853 | if (last_track==track) { |
854 | last_charge+=iqpad; |
855 | trlist->RemoveAt(last_entry); |
856 | trinfo(0)=last_track; |
857 | trinfo(1)=last_charge; |
858 | trlist->AddAt(&trinfo,last_entry); |
859 | } else { |
860 | trlist->Add(&trinfo); |
861 | } |
862 | // check the track list |
863 | Int_t nptracks=trlist->GetEntriesFast(); |
864 | if (nptracks > 2) { |
865 | for (Int_t tr=0;tr<nptracks;tr++) { |
866 | TVector *pptrk_p=(TVector*)trlist->At(tr); |
867 | TVector &pptrk=*pptrk_p; |
868 | trk[tr]=Int_t(pptrk(0)); |
869 | chtrk[tr]=Int_t(pptrk(1)); |
870 | } |
871 | } // end if nptracks |
872 | } // end if pdigit |
873 | } //end loop over clusters |
874 | } // hit loop |
875 | } // track loop |
876 | |
877 | //Int_t nentr1=list->GetEntriesFast(); |
878 | //printf(" \n counter, nentr1 %d %d\n",counter,nentr1); |
879 | |
880 | // open the file with background |
881 | |
882 | if (Add ) { |
a6f39961 |
883 | ntracks =(Int_t)TrH1->GetEntries(); |
a897a37a |
884 | //printf("background - icat,ntracks1 %d %d\n",icat,ntracks); |
885 | //printf("background - Start loop over tracks \n"); |
886 | // |
887 | // Loop over tracks |
888 | // |
889 | for (Int_t track=0; track<ntracks; track++) { |
890 | |
891 | if (fHits2) fHits2->Clear(); |
892 | if (fClusters2) fClusters2->Clear(); |
893 | |
a6f39961 |
894 | TrH1->GetEvent(track); |
a897a37a |
895 | // |
896 | // Loop over hits |
897 | AliMUONhit* mHit; |
898 | for(int i=0;i<fHits2->GetEntriesFast();++i) |
899 | { |
900 | mHit=(AliMUONhit*) (*fHits2)[i]; |
901 | Int_t nch = mHit->fChamber-1; // chamber number |
902 | if (nch >9) continue; |
903 | iChamber = &(MUON->Chamber(nch)); |
904 | Int_t rmin = (Int_t)iChamber->RInner(); |
905 | Int_t rmax = (Int_t)iChamber->ROuter(); |
906 | Float_t xbgr=mHit->fX; |
907 | Float_t ybgr=mHit->fY; |
a6f39961 |
908 | Bool_t cond=kFALSE; |
a897a37a |
909 | |
910 | for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) { |
911 | Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon]) |
912 | +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]); |
a6f39961 |
913 | if (dist<100) cond=kTRUE; |
a897a37a |
914 | } |
915 | if (!cond) continue; |
916 | |
917 | // |
918 | // Loop over pad hits |
a897a37a |
919 | for (AliMUONcluster* mPad= |
920 | (AliMUONcluster*)MUON->FirstPad(mHit,fClusters2); |
921 | mPad; |
922 | mPad=(AliMUONcluster*)MUON->NextPad(fClusters2)) |
923 | { |
e3a4d40e |
924 | |
a897a37a |
925 | Int_t cathode = mPad->fCathode; // cathode number |
926 | Int_t ipx = mPad->fPadX; // pad number on X |
927 | Int_t ipy = mPad->fPadY; // pad number on Y |
928 | Int_t iqpad = Int_t(mPad->fQpad*kScale);// charge per pad |
929 | // Int_t iqpad = mPad->fQpad; // charge per pad |
930 | |
931 | if (cathode != (icat+1)) continue; |
e3a4d40e |
932 | //if (!HitMap[nch]->CheckBoundary()) continue; |
a897a37a |
933 | // fill the info array |
934 | Float_t thex, they; |
935 | segmentation=iChamber->GetSegmentationModel(cathode); |
936 | segmentation->GetPadCxy(ipx,ipy,thex,they); |
937 | Float_t rpad=TMath::Sqrt(thex*thex+they*they); |
938 | if (rpad < rmin || iqpad ==0 || rpad > rmax) continue; |
939 | |
940 | new((*p_adr)[countadr++]) TVector(2); |
941 | TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]); |
942 | trinfo(0)=-1; // tag background |
943 | trinfo(1)=-1; |
944 | |
945 | digits[0]=ipx; |
946 | digits[1]=ipy; |
947 | digits[2]=iqpad; |
948 | digits[3]=0; |
949 | digits[4]=-1; |
950 | |
951 | AliMUONlist* pdigit; |
952 | // build the list of fired pads and update the info |
953 | if (!HitMap[nch]->TestHit(ipx, ipy)) { |
954 | list->AddAtAndExpand(new AliMUONlist(nch,digits),counter); |
955 | |
956 | HitMap[nch]->SetHit(ipx, ipy, counter); |
957 | counter++; |
958 | |
959 | pdigit=(AliMUONlist*)list->At(list->GetLast()); |
960 | // list of tracks |
961 | TObjArray *trlist=(TObjArray*)pdigit-> |
962 | TrackList(); |
963 | trlist->Add(&trinfo); |
964 | } else { |
965 | pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy); |
966 | // update charge |
967 | (*pdigit).fSignal+=iqpad; |
968 | |
969 | // update list of tracks |
970 | TObjArray* trlist=(TObjArray*)pdigit-> |
971 | TrackList(); |
972 | Int_t last_entry=trlist->GetLast(); |
973 | TVector *ptrk_p=(TVector*)trlist-> |
974 | At(last_entry); |
975 | TVector &ptrk=*ptrk_p; |
976 | Int_t last_track=Int_t(ptrk(0)); |
977 | if (last_track==-1) { |
978 | continue; |
979 | } else { |
980 | trlist->Add(&trinfo); |
981 | } |
982 | // check the track list |
983 | Int_t nptracks=trlist->GetEntriesFast(); |
984 | if (nptracks > 0) { |
985 | for (Int_t tr=0;tr<nptracks;tr++) { |
986 | TVector *pptrk_p=(TVector*)trlist->At(tr); |
987 | TVector &pptrk=*pptrk_p; |
988 | trk[tr]=Int_t(pptrk(0)); |
989 | chtrk[tr]=Int_t(pptrk(1)); |
990 | } |
991 | } // end if nptracks |
992 | } // end if pdigit |
993 | } //end loop over clusters |
994 | } // hit loop |
995 | } // track loop |
996 | //Int_t nentr2=list->GetEntriesFast(); |
997 | //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2); |
998 | TTree *fAli=gAlice->TreeK(); |
999 | TFile *file; |
1000 | |
1001 | if (fAli) file =fAli->GetCurrentFile(); |
1002 | file->cd(); |
1003 | } // if Add |
1004 | |
1005 | Int_t tracks[10]; |
1006 | Int_t charges[10]; |
1007 | //cout<<"start filling digits \n "<<endl; |
1008 | // const Float_t zero_supm = 6.; |
1009 | Int_t nentries=list->GetEntriesFast(); |
1010 | //printf(" \n \n nentries %d \n",nentries); |
1011 | // start filling the digits |
1012 | |
1013 | for (Int_t nent=0;nent<nentries;nent++) { |
1014 | AliMUONlist *address=(AliMUONlist*)list->At(nent); |
1015 | if (address==0) continue; |
1016 | Int_t ich=address->fChamber; |
1017 | Int_t q=address->fSignal; |
1018 | iChamber=(AliMUONchamber*) (*fChambers)[ich]; |
1019 | AliMUONresponse * response=iChamber->GetResponseModel(); |
1020 | Int_t adcmax= (Int_t) response->MaxAdc(); |
1021 | // add white noise and do zero-suppression and signal truncation |
1022 | Float_t MeanNoise = gRandom->Gaus(1, 0.2); |
1023 | Float_t Noise = gRandom->Gaus(0, MeanNoise); |
1024 | q+=(Int_t)Noise; |
1025 | if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise; |
1026 | if ( q <= zero_supm ) continue; |
1027 | if ( q > adcmax) q=adcmax; |
1028 | digits[0]=address->fPadX; |
1029 | digits[1]=address->fPadY; |
1030 | digits[2]=q; |
1031 | digits[3]=address->fPhysics; |
1032 | digits[4]=address->fHit; |
1033 | //printf("fSignal, fPhysics fTrack %d %d %d \n",digits[2],digits[3],digits[4]); |
1034 | |
1035 | TObjArray* trlist=(TObjArray*)address->TrackList(); |
1036 | Int_t nptracks=trlist->GetEntriesFast(); |
1037 | //printf("nptracks, trlist %d %p\n",nptracks,trlist); |
1038 | |
1039 | // this was changed to accomodate the real number of tracks |
1040 | if (nptracks > 10) { |
1041 | cout<<"Attention - nptracks > 10 "<<nptracks<<endl; |
1042 | nptracks=10; |
1043 | } |
1044 | if (nptracks > 2) { |
1045 | printf("Attention - nptracks > 2 %d \n",nptracks); |
1046 | printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q); |
1047 | } |
1048 | for (Int_t tr=0;tr<nptracks;tr++) { |
1049 | TVector *pp_p=(TVector*)trlist->At(tr); |
1050 | if(!pp_p ) printf("pp_p - %p\n",pp_p); |
1051 | TVector &pp =*pp_p; |
1052 | tracks[tr]=Int_t(pp(0)); |
1053 | charges[tr]=Int_t(pp(1)); |
1054 | //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]); |
1055 | } //end loop over list of tracks for one pad |
1056 | // Sort list of tracks according to charge |
1057 | if (nptracks > 1) { |
1058 | SortTracks(tracks,charges,nptracks); |
1059 | } |
1060 | if (nptracks < 10 ) { |
1061 | for (Int_t i=nptracks; i<10; i++) { |
1062 | tracks[i]=0; |
1063 | charges[i]=0; |
1064 | } |
1065 | } |
1066 | |
1067 | // fill digits |
1068 | MUON->AddDigits(ich,tracks,charges,digits); |
a897a37a |
1069 | } |
1070 | //cout<<"I'm out of the loops for digitisation"<<endl; |
a897a37a |
1071 | gAlice->TreeD()->Fill(); |
e3a4d40e |
1072 | TTree *TD=gAlice->TreeD(); |
1073 | |
a897a37a |
1074 | Stat_t ndig=TD->GetEntries(); |
1075 | cout<<"number of digits "<<ndig<<endl; |
1076 | TClonesArray *fDch; |
e3a4d40e |
1077 | for (int k=0;k<10;k++) { |
1078 | fDch= MUON->DigitsAddress(k); |
a897a37a |
1079 | int ndig=fDch->GetEntriesFast(); |
e3a4d40e |
1080 | printf (" i, ndig %d %d \n",k,ndig); |
a897a37a |
1081 | } |
e3a4d40e |
1082 | |
a897a37a |
1083 | MUON->ResetDigits(); |
1084 | list->Delete(); |
a6f39961 |
1085 | for(Int_t ii=0;ii<10;++ii) { |
1086 | if (HitMap[ii]) { |
1087 | hm=HitMap[ii]; |
a897a37a |
1088 | delete hm; |
a6f39961 |
1089 | HitMap[ii]=0; |
fe4da5cc |
1090 | } |
a897a37a |
1091 | } |
1092 | |
1093 | } //end loop over cathodes |
1094 | |
1095 | char hname[30]; |
1096 | sprintf(hname,"TreeD%d",nev); |
1097 | gAlice->TreeD()->Write(hname); |
1098 | // reset tree |
1099 | gAlice->TreeD()->Reset(); |
1100 | delete list; |
1101 | //Int_t nadr=p_adr->GetEntriesFast(); |
1102 | // printf(" \n \n nadr %d \n",nadr); |
1103 | |
a897a37a |
1104 | p_adr->Clear(); |
1105 | // gObjectTable->Print(); |
1106 | |
1107 | } |
1108 | |
1109 | void AliMUON::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr) |
1110 | { |
1111 | // |
1112 | // Sort the list of tracks contributing to a given digit |
1113 | // Only the 3 most significant tracks are acctually sorted |
1114 | // |
1115 | |
1116 | // |
1117 | // Loop over signals, only 3 times |
1118 | // |
1119 | |
1120 | Int_t qmax; |
1121 | Int_t jmax; |
1122 | Int_t idx[3] = {-2,-2,-2}; |
1123 | Int_t jch[3] = {-2,-2,-2}; |
1124 | Int_t jtr[3] = {-2,-2,-2}; |
1125 | Int_t i,j,imax; |
1126 | |
1127 | if (ntr<3) imax=ntr; |
1128 | else imax=3; |
1129 | for(i=0;i<imax;i++){ |
1130 | qmax=0; |
1131 | jmax=0; |
1132 | |
1133 | for(j=0;j<ntr;j++){ |
1134 | |
1135 | if((i == 1 && j == idx[i-1]) |
1136 | ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue; |
1137 | |
1138 | if(charges[j] > qmax) { |
1139 | qmax = charges[j]; |
1140 | jmax=j; |
1141 | } |
1142 | } |
1143 | |
1144 | if(qmax > 0) { |
1145 | idx[i]=jmax; |
1146 | jch[i]=charges[jmax]; |
1147 | jtr[i]=tracks[jmax]; |
1148 | } |
1149 | |
1150 | } |
1151 | |
1152 | for(i=0;i<3;i++){ |
1153 | if (jtr[i] == -2) { |
1154 | charges[i]=0; |
1155 | tracks[i]=0; |
1156 | } else { |
1157 | charges[i]=jch[i]; |
1158 | tracks[i]=jtr[i]; |
1159 | } |
1160 | } |
1161 | |
fe4da5cc |
1162 | } |
1163 | |
a897a37a |
1164 | void AliMUON::FindClusters(Int_t nev,Int_t last_entry) |
1165 | { |
1166 | |
1167 | // |
1168 | // Loop on chambers and on cathode planes |
1169 | // |
1170 | for (Int_t icat=0;icat<2;icat++) { |
1171 | gAlice->ResetDigits(); |
1172 | gAlice->TreeD()->GetEvent(last_entry+icat); // spurious +1 ... |
1173 | if (nev < 10) printf("last_entry , icat - %d %d \n",last_entry,icat); |
1174 | //gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ... |
1175 | |
1176 | for (Int_t ich=0;ich<10;ich++) { |
1177 | AliMUONchamber* iChamber=(AliMUONchamber*) (*fChambers)[ich]; |
1178 | TClonesArray *MUONdigits = this->DigitsAddress(ich); |
1179 | if (MUONdigits == 0) continue; |
1180 | // |
1181 | // Get ready the current chamber stuff |
1182 | // |
1183 | AliMUONresponse* response = iChamber->GetResponseModel(); |
1184 | AliMUONsegmentation* seg = iChamber->GetSegmentationModel(icat+1); |
1185 | AliMUONClusterFinder* rec = iChamber->GetReconstructionModel(); |
e3a4d40e |
1186 | //printf("icat, ich, seg - %d %d %p\n",icat,ich,seg); |
a897a37a |
1187 | if (seg) { |
1188 | rec->SetSegmentation(seg); |
1189 | rec->SetResponse(response); |
1190 | rec->SetDigits(MUONdigits); |
1191 | rec->SetChamber(ich); |
1192 | if (nev==0) rec->CalibrateCOG(); |
a897a37a |
1193 | rec->FindRawClusters(); |
1194 | } |
1195 | //printf("Finish FindRawClusters for cathode %d in chamber %d\n",icat,ich); |
1196 | |
1197 | TClonesArray *fRch; |
1198 | fRch=RawClustAddress(ich); |
1199 | fRch->Sort(); |
1200 | // it seems to work |
1201 | |
1202 | |
1203 | } // for ich |
1204 | // fill the tree |
e3a4d40e |
1205 | TTree *TR=gAlice->TreeR(); |
a897a37a |
1206 | |
1207 | gAlice->TreeR()->Fill(); |
e3a4d40e |
1208 | |
a897a37a |
1209 | Stat_t nent=TR->GetEntries(); |
1210 | cout<<"number of entries "<<nent<<endl; |
1211 | TClonesArray *fRch; |
1212 | for (int i=0;i<10;i++) { |
1213 | fRch=RawClustAddress(i); |
1214 | int nraw=fRch->GetEntriesFast(); |
1215 | printf (" i, nraw %d %d \n",i,nraw); |
1216 | } |
a897a37a |
1217 | ResetRawClusters(); |
1218 | |
1219 | } // for icat |
1220 | |
1221 | char hname[30]; |
1222 | sprintf(hname,"TreeR%d",nev); |
1223 | gAlice->TreeR()->Write(hname); |
1224 | gAlice->TreeR()->Reset(); |
1225 | |
1226 | //gObjectTable->Print(); |
1227 | |
1228 | } |
fe4da5cc |
1229 | |
1230 | //______________________________________________________________________________ |
a897a37a |
1231 | //_____________________________________________________________________________ |
1232 | void AliMUON::CathodeCorrelation(Int_t nev) |
1233 | { |
1234 | |
1235 | // Correlates the clusters on the two cathode planes and build a list of |
1236 | // other possible combinations (potential ghosts) - for the moment use the |
1237 | // criteria of minimum distance between the CoGs of the two correlated |
1238 | // clusters |
1239 | |
1240 | |
1241 | // |
1242 | // Loop on chambers and on clusters on the cathode plane with the highest |
1243 | // number of clusters |
1244 | |
a6f39961 |
1245 | static Bool_t first=kTRUE; |
a897a37a |
1246 | |
1247 | AliMUONRawCluster *mRaw1; |
1248 | AliMUONRawCluster *mRaw2; |
1249 | AliMUONchamber *iChamber; |
1250 | AliMUONsegmentation *seg; |
1251 | TArrayF x1, y1, x2, y2, q1, q2; |
1252 | x1.Set(5000); |
1253 | x2.Set(5000); |
1254 | y1.Set(5000); |
1255 | y2.Set(5000); |
1256 | q1.Set(5000); |
1257 | q2.Set(5000); |
1258 | |
1259 | // Get pointers to Alice detectors and Digits containers |
1260 | TTree *TR = gAlice->TreeR(); |
1261 | Int_t nent=(Int_t)TR->GetEntries(); |
1262 | if (nev < 10) printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent); |
1263 | |
1264 | |
1265 | Int_t idx[4]; |
1266 | Float_t xc2[4],yc2[4]; |
1267 | Float_t xrec2, yrec2; |
1268 | Float_t xd0, xdif, ydif; |
1269 | Float_t ysrch,xd,xmax,ymax; |
a6f39961 |
1270 | Int_t ilow, iup, iraw1, i; |
a897a37a |
1271 | // |
1272 | Float_t xarray[50]; |
1273 | Float_t xdarray[50]; |
1274 | Float_t yarray[50]; |
1275 | Float_t qarray[50]; |
1276 | Int_t idx2[50]; |
1277 | |
1278 | // Int_t nraw[2], entry,cathode; |
1279 | |
a6f39961 |
1280 | for (i=0;i<50;i++) { |
a897a37a |
1281 | xdarray[i]=1100.; |
1282 | xarray[i]=0.; |
1283 | yarray[i]=0.; |
1284 | qarray[i]=0.; |
1285 | idx2[i]=-1; |
1286 | } |
a6f39961 |
1287 | for (i=0;i<4;i++) { |
a897a37a |
1288 | idx[i]=-1; |
1289 | xc2[i]=0.; |
1290 | yc2[i]=0.; |
1291 | } |
1292 | |
1293 | // access to the Raw Clusters tree |
1294 | for (Int_t ich=0;ich<10;ich++) { |
1295 | iChamber = &(Chamber(ich)); |
1296 | TClonesArray *MUONrawclust = RawClustAddress(ich); |
1297 | ResetRawClusters(); |
1298 | TR->GetEvent(nent-2); |
1299 | //TR->GetEvent(1); |
1300 | Int_t nrawcl1 = MUONrawclust->GetEntries(); |
1301 | // printf("Found %d raw clusters for cathode 1 in chamber %d \n" |
1302 | // ,nrawcl1,ich+1); |
1303 | if (!nrawcl1) continue; |
1304 | |
1305 | seg = iChamber->GetSegmentationModel(1); |
1306 | // loop over raw clusters of first cathode |
1307 | for (iraw1=0; iraw1<nrawcl1; iraw1++) { |
1308 | mRaw1= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw1); |
1309 | x1[iraw1]=mRaw1->fX; |
1310 | y1[iraw1]=mRaw1->fY; |
1311 | q1[iraw1]=(Float_t)mRaw1->fQ; //maybe better fPeakSignal |
1312 | } // rawclusters cathode 1 |
1313 | // |
1314 | // Get information from 2nd cathode |
1315 | ResetRawClusters(); |
1316 | TR->GetEvent(nent-1); |
1317 | //TR->GetEvent(2); |
1318 | Int_t nrawcl2 = MUONrawclust->GetEntries(); |
1319 | if (!nrawcl2) { |
1320 | for (iraw1=0; iraw1<nrawcl1; iraw1++) { |
1321 | idx[3]=iraw1; |
1322 | xc2[3]=x1[iraw1]; |
1323 | yc2[3]=y1[iraw1]; |
1324 | //printf("nrawcl2 is zero - idx[0] %d\n",idx[0]); |
1325 | |
1326 | AddCathCorrel(ich,idx,xc2,yc2); |
1327 | // reset |
1328 | idx[3]=-1; |
1329 | xc2[3]=0.; |
1330 | yc2[3]=0.; |
1331 | |
1332 | } // store information from cathode 1 only |
1333 | } else { |
1334 | // printf("Found %d raw clusters for cathode 2 in chamber %d \n", |
1335 | // nrawcl2, ich+1); |
1336 | |
1337 | for (Int_t iraw2=0; iraw2<nrawcl2; iraw2++) { |
1338 | mRaw2= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw2); |
1339 | x2[iraw2]=mRaw2->fX; |
1340 | y2[iraw2]=mRaw2->fY; |
1341 | q2[iraw2]=(Float_t)mRaw2->fQ; |
1342 | } // rawclusters cathode 2 |
1343 | // |
1344 | // Initalisation finished |
1345 | for (iraw1=0; iraw1<nrawcl1; iraw1++) { |
1346 | // find the sector |
1347 | Int_t ix,iy; |
1348 | seg->GetPadIxy(x1[iraw1],y1[iraw1],ix,iy); |
1349 | Int_t isec=seg->Sector(ix,iy); |
1350 | // range to look for ghosts ?! |
1351 | if (ich < 5) { |
1352 | ymax = seg->Dpy(isec)*7/2; |
1353 | xmax = seg->Dpx(isec)*7/2; |
1354 | } else { |
1355 | ymax = seg->Dpy(isec)*13/2; |
1356 | xmax = seg->Dpx(isec)*3/2; |
1357 | } |
1358 | ysrch=ymax+y1[iraw1]; |
1359 | |
1360 | ilow = AliMUONRawCluster:: |
1361 | BinarySearch(ysrch-2*ymax,y2,0,nrawcl2+1); |
1362 | iup= AliMUONRawCluster:: |
1363 | BinarySearch(ysrch,y2,ilow,nrawcl2+1); |
1364 | if (ilow<0 || iup <0 || iup>nrawcl2) continue; |
1365 | Int_t counter=0; |
1366 | for (Int_t iraw2=ilow; iraw2<=iup; iraw2++) { |
1367 | xrec2=x2[iraw2]; |
1368 | yrec2=y2[iraw2]; |
1369 | xdif=x1[iraw1]-xrec2; |
1370 | ydif=y1[iraw1]-yrec2; |
1371 | xd=TMath::Sqrt(xdif*xdif+ydif*ydif); |
1372 | if (iraw2==ilow) { |
1373 | if (ilow==iup) |
1374 | xd0=TMath:: |
1375 | Sqrt(2*xmax*2*xmax+2*ymax*2*ymax); |
1376 | else xd0=101.; |
1377 | } |
1378 | Float_t qdif=TMath::Abs(q1[iraw1]-q2[iraw2])/q1[iraw1]; |
1379 | |
1380 | if (x1[iraw1]*xrec2 > 0) { |
1381 | if (xd <= xd0 ) { |
1382 | // printf("q1, q2 qdif % f %f %f \n",q1[iraw1],q2[iraw2],qdif); |
1383 | // printf("x1, x2 y1 y2 % f %f %f %f \n",x1[iraw1],xrec2,y1[iraw1],yrec2); |
1384 | //if (qdif <0.3) { //check this number |
1385 | |
1386 | xd0=xd; |
1387 | idx2[counter]=iraw2; |
1388 | xdarray[counter]=xd; |
1389 | xarray[counter]=xdif; |
1390 | yarray[counter]=ydif; |
1391 | qarray[counter]=qdif; |
1392 | counter++; |
1393 | // } |
1394 | |
1395 | } |
1396 | } // check for same quadrant |
1397 | } // loop over 2nd cathode range |
1398 | |
1399 | |
1400 | if (counter >=2) { |
1401 | AliMUONRawCluster:: |
1402 | SortMin(idx2,xdarray,xarray,yarray,qarray,counter); |
1403 | if (xdarray[0]<seg->Dpx(isec) && xdarray[1]<seg->Dpx(isec)) { |
1404 | if (qarray[0]>qarray[1]){ |
1405 | Int_t swap=idx2[0]; |
1406 | idx2[0]=idx2[1]; |
1407 | idx2[1]=swap; |
1408 | } |
1409 | } |
1410 | } |
1411 | int imax; |
1412 | if (counter <3) imax=counter; |
1413 | else imax=3; |
1414 | |
1415 | for (int i=0;i<imax;i++) { |
1416 | if (idx2[i] >= 0 && idx2[i] < nrawcl2) { |
1417 | if (xarray[i] > xmax || yarray[i] > 2*ymax) |
1418 | continue; |
1419 | idx[i]=idx2[i]; |
1420 | xc2[i]=x2[idx2[i]]; |
1421 | yc2[i]=y2[idx2[i]]; |
1422 | } |
1423 | } |
1424 | // add info about the cluster on the 'starting' cathode |
1425 | |
1426 | idx[3]=iraw1; |
1427 | xc2[3]=x1[iraw1]; |
1428 | yc2[3]=y1[iraw1]; |
1429 | //if (idx[0] <0) printf("iraw1 imax idx2[0] idx[0] %d %d %d %d\n",iraw1,imax,idx2[0],idx[0]); |
1430 | AddCathCorrel(ich,idx,xc2,yc2); |
1431 | // reset |
a6f39961 |
1432 | for (Int_t ii=0;ii<counter;ii++) { |
1433 | xdarray[ii]=1100.; |
1434 | xarray[ii]=0.; |
1435 | yarray[ii]=0.; |
1436 | qarray[ii]=0.; |
1437 | idx2[ii]=-1; |
a897a37a |
1438 | } |
a6f39961 |
1439 | for (Int_t iii=0;iii<3;iii++) { |
1440 | idx[iii]=-1; |
1441 | xc2[iii]=0.; |
1442 | yc2[iii]=0.; |
a897a37a |
1443 | } |
1444 | } // iraw1 |
1445 | } |
1446 | x1.Reset(); |
1447 | x2.Reset(); |
1448 | y1.Reset(); |
1449 | y2.Reset(); |
1450 | q1.Reset(); |
1451 | q2.Reset(); |
1452 | } //ich |
1453 | // |
1454 | if (first) { |
1455 | MakeTreeC("C"); |
a6f39961 |
1456 | first=kFALSE; |
a897a37a |
1457 | } |
1458 | TTree *TC=TreeC(); |
1459 | TC->Fill(); |
1460 | //Int_t nentries=(Int_t)TC->GetEntries(); |
1461 | //cout<<"number entries in tree of correlated clusters "<<nentries<<endl; |
1462 | TClonesArray *fCch; |
1463 | static Int_t countev=0; |
1464 | Int_t countch=0; |
1465 | |
a6f39961 |
1466 | for (Int_t ii=0;ii<10;ii++) { |
1467 | fCch= CathCorrelAddress(ii); |
a897a37a |
1468 | Int_t ncor=fCch->GetEntriesFast(); |
a6f39961 |
1469 | printf (" ii, ncor %d %d \n",ii,ncor); |
a897a37a |
1470 | if (ncor>=2) countch++; |
1471 | } |
1472 | |
1473 | // write |
1474 | char hname[30]; |
1475 | sprintf(hname,"TreeC%d",nev); |
1476 | TC->Write(hname); |
1477 | // reset tree |
1478 | ResetCorrelation(); |
1479 | TC->Reset(); |
1480 | |
1481 | if (countch==10) countev++; |
1482 | printf("countev - %d\n",countev); |
1483 | |
1484 | // gObjectTable->Print(); |
1485 | |
1486 | |
1487 | } |
1488 | |
1489 | |
1490 | //_____________________________________________________________________________ |
1491 | |
1492 | void AliMUON::MakeTreeC(Option_t *option) |
1493 | { |
1494 | char *C = strstr(option,"C"); |
1495 | if (C && !fTreeC) fTreeC = new TTree("TC","CathodeCorrelation"); |
1496 | |
1497 | // Create a branch for correlation |
1498 | |
1499 | const Int_t buffersize = 4000; |
1500 | char branchname[30]; |
1501 | |
1502 | // one branch for correlation per chamber |
1503 | for (int i=0; i<10 ;i++) { |
1504 | sprintf(branchname,"%sCorrelation%d",GetName(),i+1); |
1505 | |
1506 | if (fCathCorrel && fTreeC) { |
1507 | TreeC()->Branch(branchname,&((*fCathCorrel)[i]), buffersize); |
1508 | printf("Making Branch %s for correlation in chamber %d\n",branchname,i+1); |
1509 | } |
1510 | } |
1511 | } |
1512 | |
1513 | //_____________________________________________________________________________ |
1514 | void AliMUON::GetTreeC(Int_t event) |
1515 | { |
1516 | |
1517 | // set the branch address |
1518 | char treeName[20]; |
1519 | char branchname[30]; |
1520 | |
1521 | ResetCorrelation(); |
1522 | if (fTreeC) { |
1523 | delete fTreeC; |
1524 | } |
1525 | |
1526 | sprintf(treeName,"TreeC%d",event); |
1527 | fTreeC = (TTree*)gDirectory->Get(treeName); |
1528 | |
1529 | |
1530 | TBranch *branch; |
1531 | if (fTreeC) { |
1532 | for (int i=0; i<10; i++) { |
1533 | sprintf(branchname,"%sCorrelation%d",GetName(),i+1); |
1534 | if (fCathCorrel) { |
1535 | branch = fTreeC->GetBranch(branchname); |
1536 | if (branch) branch->SetAddress(&((*fCathCorrel)[i])); |
1537 | } |
1538 | } |
1539 | } else { |
1540 | printf("ERROR: cannot find CathodeCorrelation Tree for event:%d\n",event); |
1541 | } |
1542 | |
1543 | // gObjectTable->Print(); |
1544 | |
1545 | } |
1546 | |
1547 | |
fe4da5cc |
1548 | void AliMUON::Streamer(TBuffer &R__b) |
1549 | { |
1550 | // Stream an object of class AliMUON. |
1551 | AliMUONchamber *iChamber; |
1552 | AliMUONsegmentation *segmentation; |
1553 | AliMUONresponse *response; |
1554 | TClonesArray *digitsaddress; |
a897a37a |
1555 | TClonesArray *rawcladdress; |
1556 | TClonesArray *corcladdress; |
1557 | // TObjArray *clustaddress; |
fe4da5cc |
1558 | |
1559 | if (R__b.IsReading()) { |
1560 | Version_t R__v = R__b.ReadVersion(); if (R__v) { } |
1561 | AliDetector::Streamer(R__b); |
1562 | R__b >> fNclusters; |
1563 | R__b >> fClusters; // diff |
1564 | R__b >> fDchambers; |
a897a37a |
1565 | R__b >> fRawClusters; |
1566 | R__b >> fCathCorrel; |
fe4da5cc |
1567 | R__b.ReadArray(fNdch); |
a897a37a |
1568 | R__b.ReadArray(fNrawch); |
1569 | R__b.ReadArray(fNcorch); |
fe4da5cc |
1570 | // |
1571 | R__b >> fAccCut; |
1572 | R__b >> fAccMin; |
1573 | R__b >> fAccMax; |
1574 | // |
a897a37a |
1575 | // modifs perso |
1576 | R__b >> fSPxzCut; |
1577 | R__b >> fSSigmaCut; |
1578 | R__b >> fSXPrec; |
1579 | R__b >> fSYPrec; |
1580 | // |
fe4da5cc |
1581 | R__b >> fChambers; |
1582 | // Stream chamber related information |
1583 | for (Int_t i =0; i<10; i++) { |
1584 | iChamber=(AliMUONchamber*) (*fChambers)[i]; |
1585 | iChamber->Streamer(R__b); |
1586 | if (iChamber->Nsec()==1) { |
1587 | segmentation=iChamber->GetSegmentationModel(1); |
1588 | segmentation->Streamer(R__b); |
1589 | } else { |
1590 | segmentation=iChamber->GetSegmentationModel(1); |
1591 | segmentation->Streamer(R__b); |
1592 | segmentation=iChamber->GetSegmentationModel(2); |
1593 | segmentation->Streamer(R__b); |
1594 | } |
1595 | response=iChamber->GetResponseModel(); |
1596 | response->Streamer(R__b); |
1597 | digitsaddress=(TClonesArray*) (*fDchambers)[i]; |
1598 | digitsaddress->Streamer(R__b); |
a897a37a |
1599 | rawcladdress=(TClonesArray*) (*fRawClusters)[i]; |
1600 | rawcladdress->Streamer(R__b); |
1601 | corcladdress=(TClonesArray*) (*fCathCorrel)[i]; |
1602 | corcladdress->Streamer(R__b); |
fe4da5cc |
1603 | } |
1604 | |
1605 | } else { |
1606 | R__b.WriteVersion(AliMUON::IsA()); |
1607 | AliDetector::Streamer(R__b); |
1608 | R__b << fNclusters; |
1609 | R__b << fClusters; // diff |
1610 | R__b << fDchambers; |
a897a37a |
1611 | R__b << fRawClusters; |
1612 | R__b << fCathCorrel; |
fe4da5cc |
1613 | R__b.WriteArray(fNdch, 10); |
a897a37a |
1614 | R__b.WriteArray(fNrawch, 10); |
1615 | R__b.WriteArray(fNcorch, 10); |
fe4da5cc |
1616 | // |
1617 | R__b << fAccCut; |
1618 | R__b << fAccMin; |
1619 | R__b << fAccMax; |
1620 | // |
a897a37a |
1621 | // modifs perso |
1622 | R__b << fSPxzCut; |
1623 | R__b << fSSigmaCut; |
1624 | R__b << fSXPrec; |
1625 | R__b << fSYPrec; |
1626 | // |
fe4da5cc |
1627 | R__b << fChambers; |
1628 | // Stream chamber related information |
1629 | for (Int_t i =0; i<10; i++) { |
1630 | iChamber=(AliMUONchamber*) (*fChambers)[i]; |
1631 | iChamber->Streamer(R__b); |
1632 | if (iChamber->Nsec()==1) { |
1633 | segmentation=iChamber->GetSegmentationModel(1); |
1634 | segmentation->Streamer(R__b); |
1635 | } else { |
1636 | segmentation=iChamber->GetSegmentationModel(1); |
1637 | segmentation->Streamer(R__b); |
1638 | segmentation=iChamber->GetSegmentationModel(2); |
1639 | segmentation->Streamer(R__b); |
1640 | } |
1641 | response=iChamber->GetResponseModel(); |
1642 | response->Streamer(R__b); |
fe4da5cc |
1643 | digitsaddress=(TClonesArray*) (*fDchambers)[i]; |
1644 | digitsaddress->Streamer(R__b); |
a897a37a |
1645 | rawcladdress=(TClonesArray*) (*fRawClusters)[i]; |
1646 | rawcladdress->Streamer(R__b); |
1647 | corcladdress=(TClonesArray*) (*fCathCorrel)[i]; |
1648 | corcladdress->Streamer(R__b); |
fe4da5cc |
1649 | } |
1650 | } |
1651 | } |
a897a37a |
1652 | AliMUONcluster* AliMUON::FirstPad(AliMUONhit* hit, TClonesArray *clusters) |
fe4da5cc |
1653 | { |
1654 | // |
1655 | // Initialise the pad iterator |
1656 | // Return the address of the first padhit for hit |
a897a37a |
1657 | TClonesArray *theClusters = clusters; |
fe4da5cc |
1658 | Int_t nclust = theClusters->GetEntriesFast(); |
1659 | if (nclust && hit->fPHlast > 0) { |
1660 | sMaxIterPad=hit->fPHlast; |
1661 | sCurIterPad=hit->fPHfirst; |
a897a37a |
1662 | return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1); |
fe4da5cc |
1663 | } else { |
1664 | return 0; |
1665 | } |
1666 | } |
1667 | |
a897a37a |
1668 | AliMUONcluster* AliMUON::NextPad(TClonesArray *clusters) |
fe4da5cc |
1669 | { |
1670 | sCurIterPad++; |
1671 | if (sCurIterPad <= sMaxIterPad) { |
a897a37a |
1672 | return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1); |
fe4da5cc |
1673 | } else { |
1674 | return 0; |
1675 | } |
1676 | } |
1677 | |
a897a37a |
1678 | //////////////////////////// modifs perso /////////////// |
1679 | |
1680 | static TTree *ntuple_global; |
1681 | static TFile *hfile_global; |
1682 | |
1683 | // variables of the tracking ntuple |
1684 | struct { |
1685 | Int_t ievr; // number of event |
1686 | Int_t ntrackr; // number of tracks per event |
1687 | Int_t istatr[500]; // 1 = good muon, 2 = ghost, 0 = something else |
1688 | Int_t isignr[500]; // sign of the track |
1689 | Float_t pxr[500]; // x momentum of the reconstructed track |
1690 | Float_t pyr[500]; // y momentum of the reconstructed track |
1691 | Float_t pzr[500]; // z momentum of the reconstructed track |
1692 | Float_t zvr[500]; // z vertex |
1693 | Float_t chi2r[500]; // chi2 of the fit of the track with the field map |
1694 | Float_t pxv[500]; // x momentum at vertex |
1695 | Float_t pyv[500]; // y momentum at vertex |
1696 | Float_t pzv[500]; // z momentum at vertex |
1697 | } ntuple_st; |
1698 | |
1699 | AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster) |
1700 | { |
1701 | TClonesArray *MUONrawclust = RawClustAddress(ichamber); |
1702 | ResetRawClusters(); |
1703 | TTree *TR = gAlice->TreeR(); |
1704 | Int_t nent=(Int_t)TR->GetEntries(); |
1705 | TR->GetEvent(nent-2+icathod-1); |
1706 | //TR->GetEvent(icathod); |
1707 | //Int_t nrawcl = (Int_t)MUONrawclust->GetEntriesFast(); |
1708 | |
1709 | AliMUONRawCluster * mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(icluster); |
1710 | //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw); |
1711 | |
1712 | return mRaw; |
1713 | } |
1714 | |
1715 | void AliMUON::Reconst(Int_t &ifit, Int_t &idebug, Int_t bgd_ev, Int_t &nev, Int_t &idres, Int_t &ireadgeant, Option_t *option,Text_t *filename) |
1716 | { |
1717 | // |
1718 | // open kine and hits tree of background file for reconstruction of geant hits |
1719 | // call tracking fortran program |
a6f39961 |
1720 | static Bool_t first=kTRUE; |
a897a37a |
1721 | static TFile *File; |
1722 | char *Add = strstr(option,"Add"); |
1723 | |
1724 | if (Add ) { // only in case of background with geant hits |
1725 | if(first) { |
1726 | fFileName=filename; |
1727 | cout<<"filename "<<fFileName<<endl; |
1728 | File=new TFile(fFileName); |
1729 | cout<<"I have opened "<<fFileName<<" file "<<endl; |
1730 | fHits2 = new TClonesArray("AliMUONhit",1000 ); |
1731 | fParticles2 = new TClonesArray("GParticle",1000); |
a6f39961 |
1732 | first=kFALSE; |
a897a37a |
1733 | } |
1734 | File->cd(); |
1735 | if(fHits2) fHits2->Clear(); |
1736 | if(fParticles2) fParticles2->Clear(); |
a6f39961 |
1737 | if(TrH1) delete TrH1; |
1738 | TrH1=0; |
a897a37a |
1739 | if(TK1) delete TK1; |
1740 | TK1=0; |
1741 | // Get Hits Tree header from file |
1742 | char treeName[20]; |
1743 | sprintf(treeName,"TreeH%d",bgd_ev); |
a6f39961 |
1744 | TrH1 = (TTree*)gDirectory->Get(treeName); |
1745 | if (!TrH1) { |
a897a37a |
1746 | printf("ERROR: cannot find Hits Tree for event:%d\n",bgd_ev); |
1747 | } |
1748 | // set branch addresses |
1749 | TBranch *branch; |
1750 | char branchname[30]; |
1751 | sprintf(branchname,"%s",GetName()); |
a6f39961 |
1752 | if (TrH1 && fHits2) { |
1753 | branch = TrH1->GetBranch(branchname); |
a897a37a |
1754 | if (branch) branch->SetAddress(&fHits2); |
1755 | } |
a6f39961 |
1756 | TrH1->GetEntries(); |
a897a37a |
1757 | // get the Kine tree |
1758 | sprintf(treeName,"TreeK%d",bgd_ev); |
1759 | TK1 = (TTree*)gDirectory->Get(treeName); |
1760 | if (!TK1) { |
1761 | printf("ERROR: cannot find Kine Tree for event:%d\n",bgd_ev); |
1762 | } |
1763 | // set branch addresses |
1764 | if (TK1) |
1765 | TK1->SetBranchAddress("Particles", &fParticles2); |
1766 | TK1->GetEvent(0); |
1767 | |
1768 | // get back to the first file |
1769 | TTree *TK = gAlice->TreeK(); |
1770 | TFile *file1 = 0; |
1771 | if (TK) file1 = TK->GetCurrentFile(); |
1772 | file1->cd(); |
1773 | |
1774 | } // end if Add |
1775 | |
1776 | // call tracking fortran program |
1777 | reconstmuon(ifit,idebug,nev,idres,ireadgeant); |
1778 | } |
1779 | |
1780 | |
e3a4d40e |
1781 | void AliMUON::InitTracking(Double_t &seff, Double_t &sb0, Double_t &sbl3) |
a897a37a |
1782 | { |
1783 | // |
1784 | // introduce in fortran program somme parameters and cuts for tracking |
1785 | // create output file "reconst.root" (histos + ntuple) |
1786 | cutpxz(fSPxzCut); // Pxz cut (GeV/c) to begin the track finding |
1787 | sigmacut(fSSigmaCut); // Number of sigmas delimiting the searching areas |
1788 | xpreci(fSXPrec); // Chamber precision in X (cm) |
1789 | ypreci(fSYPrec); // Chamber precision in Y (cm) |
1790 | reco_init(seff,sb0,sbl3); |
1791 | } |
1792 | |
1793 | void AliMUON::FinishEvent() |
1794 | { |
1795 | TTree *TK = gAlice->TreeK(); |
1796 | TFile *file1 = 0; |
1797 | if (TK) file1 = TK->GetCurrentFile(); |
1798 | file1->cd(); |
1799 | } |
1800 | |
e3a4d40e |
1801 | void AliMUON::CloseTracking() |
a897a37a |
1802 | { |
1803 | // |
1804 | // write histos and ntuple to "reconst.root" file |
1805 | reco_term(); |
1806 | } |
1807 | |
e3a4d40e |
1808 | void chfill(Int_t &id, Float_t &x, Float_t &, Float_t &) |
a897a37a |
1809 | { |
1810 | // |
1811 | // fill histo like hfill in fortran |
1812 | char name[5]; |
1813 | sprintf(name,"h%d",id); |
1814 | TH1F *h1 = (TH1F*) gDirectory->Get(name); |
1815 | h1->Fill(x); |
1816 | } |
1817 | |
1818 | void chfill2(Int_t &id, Float_t &x, Float_t &y, Float_t &w) |
1819 | { |
1820 | // |
1821 | // fill histo like hfill2 in fortran |
1822 | char name[5]; |
1823 | sprintf(name,"h%d",id); |
1824 | TH2F *h2 = (TH2F*) gDirectory->Get(name); |
1825 | h2->Fill(x,y,w); |
1826 | } |
1827 | |
1828 | void chf1(Int_t &id, Float_t &x, Float_t &w) |
1829 | { |
1830 | // |
1831 | // fill histo like hf1 in fortran |
1832 | char name[5]; |
1833 | sprintf(name,"h%d",id); |
1834 | TH1F *h1 = (TH1F*) gDirectory->Get(name); |
1835 | h1->Fill(x,w); |
1836 | } |
1837 | |
1838 | void hist_create() |
1839 | { |
1840 | // |
1841 | // Create an output file ("reconst.root") |
1842 | // Create some histograms and an ntuple |
1843 | |
1844 | hfile_global = new TFile("reconst.root","RECREATE","Ntuple - reconstruction"); |
1845 | |
1846 | ntuple_global = new TTree("ntuple","Reconst ntuple"); |
1847 | ntuple_global->Branch("ievr",&ntuple_st.ievr,"ievr/I"); |
1848 | ntuple_global->Branch("ntrackr",&ntuple_st.ntrackr,"ntrackr/I"); |
1849 | ntuple_global->Branch("istatr",&ntuple_st.istatr[0],"istatr[500]/I"); |
1850 | ntuple_global->Branch("isignr",&ntuple_st.isignr[0],"isignr[500]/I"); |
1851 | ntuple_global->Branch("pxr",&ntuple_st.pxr[0],"pxr[500]/F"); |
1852 | ntuple_global->Branch("pyr",&ntuple_st.pyr[0],"pyr[500]/F"); |
1853 | ntuple_global->Branch("pzr",&ntuple_st.pzr[0],"pzr[500]/F"); |
1854 | ntuple_global->Branch("zvr",&ntuple_st.zvr[0],"zvr[500]/F"); |
1855 | ntuple_global->Branch("chi2r",&ntuple_st.chi2r[0],"chi2r[500]/F"); |
1856 | ntuple_global->Branch("pxv",&ntuple_st.pxv[0],"pxv[500]/F"); |
1857 | ntuple_global->Branch("pyv",&ntuple_st.pyv[0],"pyv[500]/F"); |
1858 | ntuple_global->Branch("pzv",&ntuple_st.pzv[0],"pzv[500]/F"); |
1859 | |
1860 | // test aliroot |
1861 | |
1862 | new TH1F("h100","particule id du hit geant",20,0.,20.); |
1863 | new TH1F("h101","position en x du hit geant",100,-200.,200.); |
1864 | new TH1F("h102","position en y du hit geant",100,-200.,200.); |
1865 | new TH1F("h103","chambre de tracking concernee",15,0.,14.); |
1866 | new TH1F("h104","moment ptot du hit geant",50,0.,100.); |
1867 | new TH1F("h105","px au vertex",50,0.,20.); |
1868 | new TH1F("h106","py au vertex",50,0.,20.); |
1869 | new TH1F("h107","pz au vertex",50,0.,20.); |
1870 | new TH1F("h108","position zv",50,-15.,15.); |
1871 | new TH1F("h109","position en x du hit reconstruit",100,-300.,300.); |
1872 | new TH1F("h110","position en y du hit reconstruit",100,-300.,300.); |
1873 | new TH1F("h111","delta x ",100,-0.4,0.4); |
1874 | new TH1F("h112","delta y ",100,-0.4,0.4); |
1875 | |
1876 | char hname[30]; |
1877 | char hname1[30]; |
1878 | for (int i=0;i<10;i++) { |
1879 | sprintf(hname,"deltax%d",i); |
1880 | sprintf(hname1,"h12%d",i); |
1881 | new TH1F(hname1,hname ,100,-0.4,0.4); |
1882 | sprintf(hname,"deltay%d",i); |
1883 | sprintf(hname1,"h13%d",i); |
1884 | new TH1F(hname1,hname ,100,-0.4,0.4); |
1885 | } |
1886 | new TH2F("h2000","VAR X st. 5",30,3.0,183.0,100,0.,25.); |
1887 | new TH2F("h2001","VAR Y st. 5",30,3.0,183.0,100,0.,25.); |
1888 | |
1889 | new TH2F("h2500","P vs X HHIT",30,3.0,183.0,200,0.,200.); |
1890 | new TH2F("h2501","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.); |
1891 | new TH2F("h2502","P vs X EPH2 st. 5",30,3.0,183.0,100,0.,0.000005); |
1892 | new TH2F("h2503","P vs X EAL2 st. 5",30,3.0,183.0,100,0.,0.01); |
1893 | //new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,1.5); |
1894 | new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,0.1); |
1895 | new TH2F("h2505","P vs X EYM2 st. 5",30,3.0,183.0,100,0.,30.); |
1896 | |
1897 | new TH2F("h2507","P vs X EPH st. 5",30,3.0,183.0,100,0.,0.003); |
1898 | new TH2F("h2508","P vs X EAL st. 5",30,3.0,183.0,100,0.,0.3); |
1899 | //new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,1.5); |
1900 | new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,0.4); |
1901 | new TH2F("h2510","P vs X EYM st. 5",30,3.0,183.0,100,0.,30.); |
1902 | |
1903 | new TH2F("h2511","P vs X EPH cut st. 5",30,3.0,183.0,100,0.,0.01); |
1904 | new TH2F("h2512","P vs X EAL cut st. 5",30,3.0,183.0,100,0.,0.3); |
1905 | //new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,1.5); |
1906 | new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,0.4); |
1907 | new TH2F("h2514","P vs X EYM cut st. 5",30,3.0,183.0,100,0.,30.); |
1908 | // 4 |
1909 | new TH2F("h2400","P vs X HHIT",30,3.0,183.0,200,0.,200.); |
1910 | new TH2F("h2401","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.); |
1911 | new TH2F("h2402","P vs X EPH2 st. 4",30,3.0,183.0,100,0.,0.000005); |
1912 | new TH2F("h2403","P vs X EAL2 st. 4",30,3.0,183.0,100,0.,0.05); |
1913 | //new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,1.5); |
1914 | new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,0.1); |
1915 | new TH2F("h2405","P vs X EYM2 st. 4",30,3.0,183.0,100,0.,30.); |
1916 | |
1917 | new TH2F("h2407","P vs X EPH st. 4",30,3.0,183.0,100,0.,0.003); |
1918 | new TH2F("h2408","P vs X EAL st. 4",30,3.0,183.0,100,0.,0.3); |
1919 | //new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,1.5); |
1920 | new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,0.1); |
1921 | new TH2F("h2410","P vs X EYM st. 4",30,3.0,183.0,100,0.,30.); |
1922 | |
1923 | new TH2F("h2411","P vs X EPH cut st. 4",30,3.0,183.0,100,0.,0.01); |
1924 | new TH2F("h2412","P vs X EAL cut st. 4",30,3.0,183.0,100,0.,0.3); |
1925 | //new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,1.5); |
1926 | new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,0.1); |
1927 | new TH2F("h2414","P vs X EYM cut st. 4",30,3.0,183.0,100,0.,30.); |
1928 | // 3 |
1929 | new TH1F("h2301","P2",30,3.0,183.0); |
1930 | new TH2F("h2302","P2 vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006); |
1931 | new TH2F("h2303","P2 vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.0005); |
1932 | //new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5); |
1933 | new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.); |
1934 | new TH2F("h2305","P2 vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.); |
1935 | |
1936 | new TH2F("h2307","P vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006); |
1937 | new TH2F("h2308","P vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.005); |
1938 | //new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5); |
1939 | new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.); |
1940 | new TH2F("h2310","P vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.); |
1941 | |
1942 | new TH2F("h2311","P vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06); |
1943 | new TH2F("h2312","P vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05); |
1944 | //new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5); |
1945 | new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.); |
1946 | new TH2F("h2314","P vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.); |
1947 | |
1948 | new TH2F("h2315","P2 vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06); |
1949 | new TH2F("h2316","P2 vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05); |
1950 | //new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5); |
1951 | new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.); |
1952 | new TH2F("h2318","P2 vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.); |
1953 | |
1954 | // 2 |
1955 | new TH1F("h2201","P2",30,3.0,183.0); |
1956 | new TH2F("h2202","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006); |
1957 | new TH2F("h2203","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005); |
1958 | //new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5); |
1959 | new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.); |
1960 | new TH2F("h2205","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.); |
1961 | |
1962 | new TH2F("h2207","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006); |
1963 | new TH2F("h2208","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005); |
1964 | //new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5); |
1965 | new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.); |
1966 | new TH2F("h2210","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.); |
1967 | |
1968 | new TH2F("h2211","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05); |
1969 | new TH2F("h2212","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2); |
1970 | //new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5); |
1971 | new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.); |
1972 | new TH2F("h2214","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.); |
1973 | |
1974 | new TH2F("h2215","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05); |
1975 | new TH2F("h2216","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2); |
1976 | //new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5); |
1977 | new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.); |
1978 | new TH2F("h2218","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.); |
1979 | |
1980 | // 1 |
1981 | new TH2F("h2102","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006); |
1982 | new TH2F("h2103","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005); |
1983 | //new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5); |
1984 | new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.); |
1985 | new TH2F("h2105","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.); |
1986 | |
1987 | new TH2F("h2107","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006); |
1988 | new TH2F("h2108","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005); |
1989 | //new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5); |
1990 | new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.); |
1991 | new TH2F("h2110","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.); |
1992 | |
1993 | new TH2F("h2111","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1); |
1994 | new TH2F("h2112","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2); |
1995 | //new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5); |
1996 | new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.); |
1997 | new TH2F("h2114","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.); |
1998 | |
1999 | new TH2F("h2115","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1); |
2000 | new TH2F("h2116","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2); |
2001 | //new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5); |
2002 | new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.); |
2003 | new TH2F("h2118","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.); |
2004 | |
2005 | // 2,3,4,5 |
2006 | new TH1F("h2701","P2 fit 2",30,3.0,183.0); |
2007 | new TH2F("h2702","P2 vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006); |
2008 | new TH2F("h2703","P2 vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005); |
2009 | // new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5); |
2010 | new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.); |
2011 | new TH2F("h2705","P2 vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.); |
2012 | |
2013 | new TH2F("h2707","P vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006); |
2014 | new TH2F("h2708","P vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005); |
2015 | //new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5); |
2016 | new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.); |
2017 | new TH2F("h2710","P vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.); |
2018 | |
2019 | new TH2F("h2711","P vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07); |
2020 | new TH2F("h2712","P vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2); |
2021 | //new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5); |
2022 | new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.); |
2023 | new TH2F("h2714","P vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.); |
2024 | |
2025 | new TH2F("h2715","P2 vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07); |
2026 | new TH2F("h2716","P2 vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2); |
2027 | //new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5); |
2028 | new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.); |
2029 | new TH2F("h2718","P2 vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.); |
2030 | |
2031 | // 1,3,4,5 |
2032 | new TH1F("h2801","P2 fit 1",30,3.0,183.0); |
2033 | new TH2F("h2802","P2 vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006); |
2034 | new TH2F("h2803","P2 vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005); |
2035 | //new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5); |
2036 | new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.); |
2037 | new TH2F("h2805","P2 vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.); |
2038 | |
2039 | new TH2F("h2807","P vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006); |
2040 | new TH2F("h2808","P vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005); |
2041 | //new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5); |
2042 | new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.); |
2043 | new TH2F("h2810","P vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.); |
2044 | |
2045 | new TH2F("h2811","P vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05); |
2046 | new TH2F("h2812","P vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2); |
2047 | //new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5); |
2048 | new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.); |
2049 | new TH2F("h2814","P vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.); |
2050 | |
2051 | new TH2F("h2815","P2 vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05); |
2052 | new TH2F("h2816","P2 vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2); |
2053 | //new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5); |
2054 | new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.); |
2055 | new TH2F("h2818","P2 vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.); |
2056 | // fin de test |
2057 | |
2058 | new TH1F("h500","Acceptance en H st. 4",500,0.,500.); |
2059 | new TH1F("h600","Acceptance en H st. 5",500,0.,500.); |
2060 | new TH1F("h700","X vertex track found",200,-10.,10.); |
2061 | new TH1F("h701","Y vertex track found",200,-10.,10.); |
2062 | new TH1F("h800","Rap. muon gen.",100,0.,5.); |
2063 | new TH1F("h801","Rap. muon gen. recons.",100,0.,5.); |
2064 | new TH1F("h802","Rap. muon gen. ghost ",100,0.,5.); |
2065 | new TH1F("h900","Pt muon gen.",100,0.,20.); |
2066 | new TH1F("h901","Pt muon gen. recons.",100,0.,20.); |
2067 | new TH1F("h902","Pt muon gen. ghost",100,0.,20.); |
2068 | new TH1F("h910","phi muon gen.",100,-10.,10.); |
2069 | new TH1F("h911","phi muon gen. recons.",100,-10.,10.); |
2070 | new TH1F("h912","phi muon gen. ghost",100,-10.,10.); |
2071 | new TH2F("h1001","Y VS X hit st. 1",300,-300.,300.,300,-300.,300.); |
2072 | new TH2F("h1002","Y VS X hit st. 2",300,-300.,300.,300,-300.,300.); |
2073 | new TH2F("h1003","Y VS X hit st. 3",300,-300.,300.,300,-300.,300.); |
2074 | new TH2F("h1004","Y VS X hit st. 4",300,-300.,300.,300,-300.,300.); |
2075 | new TH2F("h1005","Y VS X hit st. 5",300,-300.,300.,300,-300.,300.); |
2076 | // Histos variance dans 4 |
2077 | new TH2F("h11","VAR X st. 4",30,3.0,183.0,100,0.,2.); |
2078 | new TH2F("h12","VAR Y st. 4",30,3.0,183.0,100,0.,600.); |
2079 | new TH2F("h13","VAR PHI st. 4",30,3.0,183.0,100,0.,0.0001); |
2080 | new TH2F("h14","VAR ALM st. 4",30,3.0,183.0,100,0.,0.05); |
2081 | new TH1F("h15","P",30,3.0,183.0); |
2082 | new TH1F("h411","VAR X st. 4",100,-1.42,1.42); |
2083 | new TH1F("h412","VAR Y st. 4",100,-25.,25.); |
2084 | new TH1F("h413","VAR PHI st. 4",100,-0.01,0.01); |
2085 | new TH1F("h414","VAR ALM st. 4",100,-0.23,0.23); |
2086 | // histo2 |
2087 | new TH2F("h211","histo2-VAR X st. 4",30,3.0,183.0,100,0.,2.); |
2088 | new TH2F("h212","histo2-VAR Y st. 4",30,3.0,183.0,100,0.,600.); |
2089 | new TH1F("h213","histo2-VAR X st. 4",100,-1.42,1.42); |
2090 | new TH1F("h214","histo2-VAR Y st. 4",100,-25.,25.); |
2091 | new TH1F("h215","histo2-P",30,3.0,183.0); |
2092 | |
2093 | // Histos variance dans 2 |
2094 | new TH2F("h21","VAR X st. 2",30,3.0,183.0,100,0.,3.); |
2095 | new TH2F("h22","VAR Y st. 2",30,3.0,183.0,100,0.,7.); |
2096 | new TH2F("h23","VAR PHI st. 2",30,3.0,183.0,100,0.,0.006); |
2097 | new TH2F("h24","VAR ALM st. 2",30,3.0,183.0,100,0.,0.005); |
2098 | new TH1F("h25","P",30,3.0,183.0); |
2099 | new TH1F("h421","VAR X st. 2",100,-1.72,1.72); |
2100 | new TH1F("h422","VAR Y st. 2",100,-2.7,2.7); |
2101 | new TH1F("h423","VAR PHI st. 2",100,-0.08,0.08); |
2102 | new TH1F("h424","VAR ALM st. 2",100,-0.072,0.072); |
2103 | // histo2 |
2104 | new TH2F("h221","histo2-VAR X st. 2",30,3.0,183.0,100,0.,3.); |
2105 | new TH2F("h222","histo2-VAR Y st. 2",30,3.0,183.0,100,0.,7.); |
2106 | new TH1F("h223","histo2-VAR X st. 2",100,-1.72,1.72); |
2107 | new TH1F("h224","histo2-VAR Y st. 2",100,-2.7,2.7); |
2108 | new TH1F("h225","histo2-P",30,3.0,183.0); |
2109 | |
2110 | // Histos variance dans 1 |
2111 | new TH2F("h31","VAR X st. 1",30,3.0,183.0,100,0.,2.); |
2112 | new TH2F("h32","VAR Y st. 1",30,3.0,183.0,100,0.,0.5); |
2113 | new TH2F("h33","VAR PHI st. 1",30,3.0,183.0,100,0.,0.006); |
2114 | new TH2F("h34","VAR ALM st. 1",30,3.0,183.0,100,0.,0.005); |
2115 | new TH1F("h35","P",30,3.0,183.0); |
2116 | new TH1F("h431","VAR X st. 1",100,-1.42,1.42); |
2117 | new TH1F("h432","VAR Y st. 1",100,-0.72,0.72); |
2118 | new TH1F("h433","VAR PHI st. 1",100,-0.08,0.08); |
2119 | new TH1F("h434","VAR ALM st. 1",100,-0.072,0.072); |
2120 | // Histos variance dans 1 |
2121 | new TH2F("h41","VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.); |
2122 | new TH2F("h42","VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.); |
2123 | new TH2F("h43","VAR PHI st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005); |
2124 | new TH2F("h44","VAR ALM st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005); |
2125 | new TH1F("h45","P",30,3.0,183.0); |
2126 | new TH1F("h441","VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.); |
2127 | new TH1F("h442","VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5); |
2128 | new TH1F("h443","VAR PHI st. 1 fit 5,4,3,2,V",100,-0.072,0.072); |
2129 | new TH1F("h444","VAR ALM st. 1 fit 5,4,3,2,V",100,-0.072,0.072); |
2130 | // histo2 |
2131 | new TH2F("h241","histo2-VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.); |
2132 | new TH2F("h242","histo2-VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.); |
2133 | new TH1F("h243","histo2-VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.); |
2134 | new TH1F("h244","histo2-VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5); |
2135 | new TH1F("h245","histo2-P",30,3.0,183.0); |
2136 | |
2137 | // Histos variance dans 2 |
2138 | new TH2F("h51","VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5); |
2139 | new TH2F("h52","VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.); |
2140 | new TH2F("h53","VAR PHI st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.005); |
2141 | new TH2F("h54","VAR ALM st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.01); |
2142 | new TH1F("h55","P",30,3.0,183.0); |
2143 | new TH1F("h451","VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72); |
2144 | new TH1F("h452","VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42); |
2145 | new TH1F("h453","VAR PHI st. 2 fit 5,4,3,1,V",100,-0.072,0.072); |
2146 | new TH1F("h454","VAR ALM st. 2 fit 5,4,3,1,V",100,-0.1,0.1); |
2147 | new TH1F("h999","PTOT",30,3.0,183.0); |
2148 | // histo2 |
2149 | new TH2F("h251","histo2-VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5); |
2150 | new TH2F("h252","histo2-VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.); |
2151 | new TH1F("h253","histo2-VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72); |
2152 | new TH1F("h254","histo2-VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42); |
2153 | new TH1F("h255","histo2-P",30,3.0,183.0); |
2154 | // Histos variance dans 3 |
2155 | new TH2F("h61","VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.); |
2156 | new TH2F("h62","VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.); |
2157 | new TH2F("h63","VAR PHI st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006); |
2158 | new TH2F("h64","VAR ALM st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006); |
2159 | new TH1F("h65","P",30,3.0,183.0); |
2160 | new TH1F("h461","VAR X st. 3 fit 4,5,V",100,-2.25,2.25); |
2161 | new TH1F("h462","VAR Y st. 3 fit 4,5,V",100,-1.42,1.42); |
2162 | new TH1F("h463","VAR PHI st. 3 fit 4,5,V",100,-0.024,0.024); |
2163 | new TH1F("h464","VAR ALM st. 3 fit 4,5,V",100,-0.024,0.024); |
2164 | // histo2 |
2165 | new TH2F("h261","histo2-VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.); |
2166 | new TH2F("h262","histo2-VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.); |
2167 | new TH1F("h263","histo2-VAR X st. 3 fit 4,5,V",100,-2.25,2.25); |
2168 | new TH1F("h264","histo2-VAR Y st. 3 fit 4,5,V",100,-1.42,1.42); |
2169 | new TH1F("h265","Phisto2-",30,3.0,183.0); |
2170 | // Histos dx,dy distribution between chambers inside stations |
2171 | new TH1F("h71","DX in st. ID-70",100,-5.,5.); |
2172 | new TH1F("h81","DY in st. ID-80",100,-5.,5.); |
2173 | new TH1F("h72","DX in st. ID-70",100,-5.,5.); |
2174 | new TH1F("h82","DY in st. ID-80",100,-5.,5.); |
2175 | new TH1F("h73","DX in st. ID-70",100,-5.,5.); |
2176 | new TH1F("h83","DY in st. ID-80",100,-5.,5.); |
2177 | new TH1F("h74","DX in st. ID-70",100,-5.,5.); |
2178 | new TH1F("h84","DY in st. ID-80",100,-5.,5.); |
2179 | new TH1F("h75","DX in st. ID-70",100,-5.,5.); |
2180 | new TH1F("h85","DY in st. ID-80",100,-5.,5.); |
2181 | } |
2182 | |
2183 | void chfnt(Int_t &ievr, Int_t &ntrackr, Int_t *istatr, Int_t *isignr, Float_t *pxr, Float_t *pyr, Float_t *pzr, Float_t *zvr, Float_t *chi2r, Float_t *pxv, Float_t *pyv, Float_t *pzv) |
2184 | { |
2185 | // |
2186 | // fill the ntuple |
2187 | ntuple_st.ievr = ievr; |
2188 | ntuple_st.ntrackr = ntrackr; |
2189 | for (Int_t i=0; i<500; i++) { |
2190 | ntuple_st.istatr[i] = istatr[i]; |
2191 | ntuple_st.isignr[i] = isignr[i]; |
2192 | ntuple_st.pxr[i] = pxr[i]; |
2193 | ntuple_st.pyr[i] = pyr[i]; |
2194 | ntuple_st.pzr[i] = pzr[i]; |
2195 | ntuple_st.zvr[i] = zvr[i]; |
2196 | ntuple_st.chi2r[i] = chi2r[i]; |
2197 | ntuple_st.pxv[i] = pxv[i]; |
2198 | ntuple_st.pyv[i] = pyv[i]; |
2199 | ntuple_st.pzv[i] = pzv[i]; |
2200 | } |
2201 | ntuple_global->Fill(); |
2202 | } |
2203 | |
2204 | void hist_closed() |
2205 | { |
2206 | // |
2207 | // write histos and ntuple to "reconst.root" file |
2208 | hfile_global->Write(); |
2209 | } |
2210 | |
2211 | void trackf_read_geant(Int_t *itypg, Double_t *xtrg, Double_t *ytrg, Double_t *ptotg, Int_t *idg, Int_t *izch, Double_t *pvert1g, Double_t *pvert2g, Double_t *pvert3g, Double_t *zvertg, Int_t &nhittot1, Double_t *cx, Double_t *cy, Double_t *cz, Int_t &ievr,Int_t &nev,Double_t *xgeant, Double_t *ygeant,Double_t *clsize1, Double_t *clsize2) |
2212 | { |
2213 | // |
2214 | // introduce aliroot variables in fortran common |
2215 | // tracking study from geant hits |
2216 | // |
2217 | |
2218 | AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); |
2219 | |
2220 | // TTree *TK = gAlice->TreeK(); |
2221 | TTree *TH = gAlice->TreeH(); |
2222 | Int_t ntracks = (Int_t)TH->GetEntries(); |
2223 | cout<<"ntrack="<<ntracks<<endl; |
2224 | |
2225 | Int_t maxidg = 0; |
2226 | Int_t nres=0; |
2227 | |
2228 | // |
2229 | // Loop over tracks |
2230 | // |
2231 | |
2232 | for (Int_t track=0; track<ntracks;track++) { |
2233 | gAlice->ResetHits(); |
2234 | TH->GetEvent(track); |
2235 | |
2236 | if (MUON) { |
2237 | // |
2238 | // Loop over hits |
2239 | // |
2240 | for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); |
2241 | mHit; |
2242 | mHit=(AliMUONhit*)MUON->NextHit()) |
2243 | { |
2244 | if (maxidg<=20000) { |
2245 | |
2246 | if (mHit->fChamber > 10) continue; |
2247 | TClonesArray *fPartArray = gAlice->Particles(); |
2248 | TParticle *Part; |
2249 | Int_t ftrack = mHit->fTrack; |
2250 | Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode(); |
2251 | |
2252 | if (id==kMuonPlus||id==kMuonMinus) { |
2253 | |
2254 | // inversion de x et y car le champ est inverse dans le programme tracking |
2255 | xtrg[maxidg] = 0; |
2256 | ytrg[maxidg] = 0; |
2257 | xgeant[maxidg] = mHit->fY; // x-pos of hit |
2258 | ygeant[maxidg] = mHit->fX; // y-pos of hit |
2259 | clsize1[maxidg] = 0; // cluster size on 1-st cathode |
2260 | clsize2[maxidg] = 0; // cluster size on 2-nd cathode |
2261 | cx[maxidg] = mHit->fCyHit; // Px/P of hit |
2262 | cy[maxidg] = mHit->fCxHit; // Py/P of hit |
2263 | cz[maxidg] = mHit->fCzHit; // Pz/P of hit |
2264 | izch[maxidg] = mHit->fChamber; |
2265 | /* |
2266 | Int_t pdgtype = Int_t(mHit->fParticle); // particle number |
2267 | itypg[maxidg] = gMC->IdFromPDG(pdgtype); |
2268 | |
2269 | */ |
2270 | if (id==kMuonPlus) itypg[maxidg] = 5; |
2271 | else itypg[maxidg] = 6; |
2272 | |
a897a37a |
2273 | ptotg[maxidg] = mHit->fPTot; // P of hit |
2274 | |
2275 | Part = (TParticle*) fPartArray->UncheckedAt(ftrack); |
2276 | Float_t thet = Part->Theta(); |
2277 | thet = thet*180./3.1416; |
2278 | |
a897a37a |
2279 | Int_t iparent = Part->GetFirstMother(); |
2280 | if (iparent >= 0) { |
2281 | Int_t ip; |
2282 | while(1) { |
2283 | ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother(); |
2284 | if (ip < 0) { |
2285 | break; |
2286 | } else { |
2287 | iparent = ip; |
2288 | } |
2289 | } |
2290 | } |
2291 | //printf("iparent - %d\n",iparent); |
2292 | Int_t id1 = ftrack; // numero de la particule generee au vertex |
2293 | Int_t idum = track+1; |
2294 | Int_t id2 = ((TParticle*) fPartArray->UncheckedAt(iparent))->GetPdgCode(); |
2295 | |
2296 | if (id2==443) id2=114; |
2297 | else id2=116; |
2298 | |
2299 | if (id2==116) { |
2300 | nres++; |
2301 | } |
2302 | //printf("id2 %d\n",id2); |
2303 | idg[maxidg] = 30000*id1+10000*idum+id2; |
2304 | |
2305 | pvert1g[maxidg] = Part->Py(); // Px vertex |
2306 | pvert2g[maxidg] = Part->Px(); // Py vertex |
2307 | pvert3g[maxidg] = Part->Pz(); // Pz vertex |
2308 | zvertg[maxidg] = Part->Vz(); // z vertex |
a897a37a |
2309 | maxidg ++; |
2310 | |
2311 | } |
2312 | } |
2313 | } // hit loop |
2314 | } // if MUON |
2315 | } // track loop first file |
2316 | |
a6f39961 |
2317 | if (TrH1 && fHits2 ) { // if background file |
2318 | ntracks =(Int_t)TrH1->GetEntries(); |
a897a37a |
2319 | printf("Trackf_read - 2-nd file - ntracks %d\n",ntracks); |
2320 | |
2321 | // Loop over tracks |
2322 | for (Int_t track=0; track<ntracks; track++) { |
2323 | |
2324 | if (fHits2) fHits2->Clear(); |
a6f39961 |
2325 | TrH1->GetEvent(track); |
a897a37a |
2326 | |
2327 | // Loop over hits |
2328 | for (int i=0;i<fHits2->GetEntriesFast();i++) |
2329 | { |
2330 | AliMUONhit *mHit=(AliMUONhit*) (*fHits2)[i]; |
2331 | |
2332 | if (mHit->fChamber > 10) continue; |
2333 | |
2334 | if (maxidg<=20000) { |
2335 | |
2336 | // inversion de x et y car le champ est inverse dans le programme tracking !!!! |
2337 | xtrg[maxidg] = 0; // only for reconstructed point |
2338 | ytrg[maxidg] = 0; // only for reconstructed point |
2339 | xgeant[maxidg] = mHit->fY; // x-pos of hit |
2340 | ygeant[maxidg] = mHit->fX; // y-pos of hit |
2341 | clsize1[maxidg] = 0; // cluster size on 1-st cathode |
2342 | clsize2[maxidg] = 0; // cluster size on 2-nd cathode |
2343 | cx[maxidg] = mHit->fCyHit; // Px/P of hit |
2344 | cy[maxidg] = mHit->fCxHit; // Py/P of hit |
2345 | cz[maxidg] = mHit->fCzHit; // Pz/P of hit |
2346 | izch[maxidg] = mHit->fChamber; // chamber number |
2347 | ptotg[maxidg] = mHit->fPTot; // P of hit |
2348 | |
2349 | Int_t ftrack = mHit->fTrack; |
2350 | Int_t id1 = ftrack; // track number |
2351 | Int_t idum = track+1; |
2352 | |
2353 | TClonesArray *fPartArray = fParticles2; |
2354 | TParticle *Part; |
e3a4d40e |
2355 | Part = (TParticle*) fPartArray->UncheckedAt(ftrack); |
a897a37a |
2356 | Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode(); |
2357 | if (id==kMuonPlus||id==kMuonMinus) { |
2358 | if (id==kMuonPlus) itypg[maxidg] = 5; |
2359 | else itypg[maxidg] = 6; |
2360 | } else itypg[maxidg]=0; |
2361 | |
2362 | Int_t id2=0; // set parent to 0 for background !! |
2363 | idg[maxidg] = 30000*id1+10000*idum+id2; |
2364 | |
2365 | pvert1g[maxidg] = Part->Py(); // Px vertex |
2366 | pvert2g[maxidg] = Part->Px(); // Py vertex |
2367 | pvert3g[maxidg] = Part->Pz(); // Pz vertex |
2368 | zvertg[maxidg] = Part->Vz(); // z vertex |
2369 | |
2370 | maxidg ++; |
2371 | |
2372 | } // check limits (maxidg) |
2373 | } // hit loop |
2374 | } // track loop |
a6f39961 |
2375 | } // if TrH1 |
a897a37a |
2376 | |
2377 | ievr = nev; |
2378 | nhittot1 = maxidg ; |
2379 | cout<<"nhittot1="<<nhittot1<<endl; |
2380 | |
2381 | static Int_t nbres=0; |
2382 | if (nres>=19) nbres++; |
2383 | printf("nres, nbres %d %d \n",nres,nbres); |
2384 | |
2385 | hfile_global->cd(); |
2386 | |
2387 | } |
2388 | |
2389 | |
2390 | |
2391 | void trackf_read_spoint(Int_t *itypg, Double_t *xtrg, Double_t *ytrg, Double_t *ptotg, Int_t *idg, Int_t *izch, Double_t *pvert1g, Double_t *pvert2g, Double_t *pvert3g, Double_t *zvertg, Int_t &nhittot1, Double_t *cx, Double_t *cy, Double_t *cz, Int_t &ievr,Int_t &nev,Double_t *xgeant, Double_t *ygeant,Double_t *clsize1, Double_t *clsize2) |
2392 | |
2393 | { |
2394 | // |
2395 | // introduce aliroot variables in fortran common |
2396 | // tracking study from reconstructed points |
2397 | // |
2398 | AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); |
2399 | |
2400 | cout<<"numero de l'evenement "<<nev<<endl; |
2401 | |
2402 | MUON->GetTreeC(nev); |
2403 | TTree *TC=MUON->TreeC(); |
2404 | TC->GetEntries(); |
2405 | |
2406 | Int_t maxidg = 0; |
2407 | Int_t nres=0; |
2408 | Int_t nncor=0; |
2409 | static Int_t nuncor=0; |
2410 | static Int_t nbadcor=0; |
2411 | AliMUONRawCluster * mRaw; |
2412 | AliMUONRawCluster * mRaw1; |
2413 | TTree *TH = gAlice->TreeH(); |
2414 | |
2415 | Int_t ihit; |
2416 | Int_t mult1, mult2; |
2417 | if (MUON) { |
2418 | for (Int_t ich=0;ich<10;ich++) { |
2419 | TClonesArray *MUONcorrel = MUON->CathCorrelAddress(ich); |
2420 | MUON->ResetCorrelation(); |
2421 | TC->GetEvent(); |
2422 | Int_t ncor = (Int_t)MUONcorrel->GetEntries(); |
2423 | if (ncor>=2) nncor++; |
2424 | if (!ncor) continue; |
2425 | |
2426 | // Loop over correlated clusters |
2427 | for (Int_t icor=0;icor<ncor;icor++) { |
2428 | AliMUONcorrelation * mCor = (AliMUONcorrelation*)MUONcorrel->UncheckedAt(icor); |
2429 | |
2430 | Int_t flag=0; // = 1 if no information in the second cathode |
2431 | Int_t index = mCor->fCorrelIndex[0]; // for the second cathode |
2432 | if (index >= 0) { |
2433 | Int_t index1 = mCor->fCorrelIndex[3]; // for the 1-st cathode |
2434 | mRaw1 = MUON->RawCluster(ich,1,index1); |
2435 | mult1=mRaw1->fMultiplicity; |
2436 | mRaw = MUON->RawCluster(ich,2,index); |
2437 | mult2=mRaw->fMultiplicity; |
2438 | } else { |
2439 | index = mCor->fCorrelIndex[3]; |
2440 | mRaw = MUON->RawCluster(ich,1,index); |
2441 | mult1=mRaw->fMultiplicity; |
2442 | mult2=0; |
2443 | flag=1; |
2444 | nuncor++; |
2445 | } |
2446 | if (!mRaw) continue; |
2447 | |
2448 | Int_t ftrack1 = mRaw->fTracks[1]; // qui doit etre le meme pour |
2449 | // la cathode 1 et 2 |
2450 | ihit= mRaw->fTracks[0]; |
2451 | //printf("icor, ftrack1 ihit %d %d %d\n",icor,ftrack1,ihit); |
2452 | |
2453 | if (mRaw->fClusterType == 0 ) { |
2454 | |
2455 | if (maxidg<=20000) { |
2456 | if (flag == 0) { |
2457 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2458 | ytrg[maxidg] = (Double_t) mCor->fX[0]; |
2459 | Int_t index1 = mCor->fCorrelIndex[3]; |
2460 | mRaw1 = MUON->RawCluster(ich,1,index1); |
2461 | if (mRaw1->fClusterType==1 || mRaw1->fClusterType==2) { |
2462 | Float_t xclust=mCor->fX[3]; |
2463 | Float_t yclust=mCor->fY[3]; |
2464 | AliMUONchamber *iChamber=&(MUON->Chamber(ich)); |
2465 | AliMUONsegmentation *seg = iChamber->GetSegmentationModel(1); |
2466 | Int_t ix,iy; |
2467 | seg->GetPadIxy(xclust,yclust,ix,iy); |
2468 | Int_t isec=seg->Sector(ix,iy); |
2469 | printf("nev, CORRELATION with pure background in chamber sector %d %d %d !!!!!!!!!!!!!!!!!!!!!\n",nev,ich+1,isec); |
2470 | nbadcor++; |
2471 | |
2472 | } // end if cluster type on cathode 1 |
2473 | }else { |
2474 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2475 | ytrg[maxidg] = (Double_t) mCor->fX[3]; |
2476 | } // if iflag |
2477 | izch[maxidg] = ich+1; |
2478 | xgeant[maxidg] = 0; |
2479 | ygeant[maxidg] = 0; |
2480 | clsize1[maxidg] = mult1; |
2481 | clsize2[maxidg] = mult2; |
2482 | |
2483 | cx[maxidg] = 0; // Px/P of hit |
2484 | cy[maxidg] = 0; // Py/P of hit |
2485 | cz[maxidg] = 0; // Pz/P of hit |
2486 | itypg[maxidg] = 0; // particle number |
2487 | ptotg[maxidg] = 0; // P of hit |
2488 | idg[maxidg] = 0; |
2489 | pvert1g[maxidg] = 0; // Px vertex |
2490 | pvert2g[maxidg] = 0; // Py vertex |
2491 | pvert3g[maxidg] = 0; // Pz vertex |
2492 | zvertg[maxidg] = 0; // z vertex |
2493 | maxidg++; |
2494 | |
2495 | }// fin maxidg |
2496 | |
2497 | } else if (mRaw->fClusterType ==1 && ftrack1 < 0) // background + resonance |
2498 | { |
2499 | nres++; |
2500 | // get indexmap and loop over digits to find the signal |
2501 | Int_t nent=(Int_t)gAlice->TreeD()->GetEntries(); |
2502 | gAlice->ResetDigits(); |
2503 | if (flag==0) { |
2504 | //gAlice->TreeD()->GetEvent(2); // cathode 2 |
2505 | gAlice->TreeD()->GetEvent(nent-1); // cathode 2 |
2506 | } else { |
2507 | //gAlice->TreeD()->GetEvent(1); // cathode 1 |
2508 | gAlice->TreeD()->GetEvent(nent-2); // cathode 1 |
2509 | } |
2510 | |
2511 | TClonesArray *MUONdigits = MUON->DigitsAddress(ich); |
2512 | Int_t mul=mRaw->fMultiplicity; |
2513 | Int_t trsign; |
2514 | for (int i=0;i<mul;i++) { |
2515 | Int_t idx=mRaw->fIndexMap[i]; |
2516 | AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx); |
2517 | trsign=dig->fTracks[0]; |
2518 | ihit=dig->fHit-1; |
2519 | if (trsign > 0 && ihit >= 0) break; |
2520 | |
2521 | } // loop over indexmap |
2522 | |
2523 | //printf("trsign, ihit %d %d\n",trsign,ihit); |
2524 | //printf("signal+background : trsign %d\n",trsign); |
2525 | |
2526 | if (trsign < 0 || ihit < 0) { // no signal muon was found |
2527 | |
2528 | if (maxidg<=20000) { |
2529 | if (flag == 0) { |
2530 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2531 | ytrg[maxidg] = (Double_t) mCor->fX[0]; |
2532 | }else { |
2533 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2534 | ytrg[maxidg] = (Double_t) mCor->fX[3]; |
2535 | } |
2536 | |
2537 | izch[maxidg] = ich+1; |
2538 | |
2539 | // initialisation of informations which |
2540 | // can't be reached for background |
2541 | |
2542 | xgeant[maxidg] = 0; // only for resonances |
2543 | ygeant[maxidg] = 0; // only for resonances |
2544 | clsize1[maxidg] = mult1; |
2545 | clsize2[maxidg] = mult2; |
2546 | |
2547 | cx[maxidg] = 0; // Px/P of hit |
2548 | cy[maxidg] = 0; // Py/P of hit |
2549 | cz[maxidg] = 0; // Pz/P of hit |
2550 | itypg[maxidg] = 0; // particle number |
2551 | ptotg[maxidg] = 0; // P of hit |
2552 | idg[maxidg] = 0; |
2553 | pvert1g[maxidg] = 0; // Px vertex |
2554 | pvert2g[maxidg] = 0; // Py vertex |
2555 | pvert3g[maxidg] = 0; // Pz vertex |
2556 | zvertg[maxidg] = 0; |
2557 | maxidg++; |
2558 | |
2559 | }// fin maxidg |
2560 | } else { // signal muon - retrieve info |
2561 | //printf("inside trsign, ihit %d %d\n",trsign,ihit); |
2562 | if (maxidg<=20000) { |
2563 | if (flag == 0) { |
2564 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2565 | ytrg[maxidg] = (Double_t) mCor->fX[0]; |
2566 | }else { |
2567 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2568 | ytrg[maxidg] = (Double_t) mCor->fX[3]; |
2569 | } |
2570 | izch[maxidg] = ich+1; |
2571 | clsize1[maxidg] = mult1; |
2572 | clsize2[maxidg] = mult2; |
2573 | |
2574 | // initialise and set to the correct values |
2575 | // if signal muons |
2576 | |
2577 | xgeant[maxidg] = 0; // only for resonances |
2578 | ygeant[maxidg] = 0; // only for resonances |
2579 | |
2580 | cx[maxidg] = 0; // Px/P of hit |
2581 | cy[maxidg] = 0; // Py/P of hit |
2582 | cz[maxidg] = 0; // Pz/P of hit |
2583 | itypg[maxidg] = 0; // particle number |
2584 | ptotg[maxidg] = 0; // P of hit |
2585 | idg[maxidg] = 0; |
2586 | pvert1g[maxidg] = 0; // Px vertex |
2587 | pvert2g[maxidg] = 0; // Py vertex |
2588 | pvert3g[maxidg] = 0; // Pz vertex |
2589 | zvertg[maxidg] = 0; |
2590 | // try to retrieve info about signal muons |
2591 | gAlice->ResetHits(); |
2592 | TH->GetEvent(trsign); |
2593 | |
2594 | TClonesArray *MUONhits = MUON->Hits(); |
2595 | AliMUONhit *mHit= (AliMUONhit*)MUONhits-> |
2596 | UncheckedAt(ihit); |
2597 | TClonesArray *fPartArray = gAlice->Particles(); |
2598 | TParticle *Part; |
2599 | Int_t nch=mHit->fChamber-1; |
2600 | //printf("sig+bgr ich, nch %d %d \n",ich,nch); |
2601 | if (nch==ich) { |
2602 | Int_t ftrack = mHit->fTrack; |
2603 | Int_t id = ((TParticle*) fPartArray-> |
2604 | UncheckedAt(ftrack))->GetPdgCode(); |
2605 | if (id==kMuonPlus||id==kMuonMinus) { |
2606 | xgeant[maxidg] = (Double_t) mHit->fY; |
2607 | ygeant[maxidg] = (Double_t) mHit->fX; |
2608 | cx[maxidg] = (Double_t) mHit->fCyHit; |
2609 | cy[maxidg] = (Double_t) mHit->fCxHit; |
2610 | cz[maxidg] = (Double_t) mHit->fCzHit; |
2611 | |
2612 | if (id==kMuonPlus) { |
2613 | itypg[maxidg] = 5; |
2614 | } else if (id==kMuonMinus) { |
2615 | itypg[maxidg] = 6; |
2616 | } else itypg[maxidg] = 0; |
2617 | |
2618 | ptotg[maxidg] = (Double_t) mHit->fPTot; |
2619 | Part = (TParticle*) fPartArray-> |
2620 | UncheckedAt(ftrack); |
2621 | Int_t iparent = Part->GetFirstMother(); |
2622 | Int_t id2; |
2623 | id2 = ((TParticle*) fPartArray-> |
2624 | UncheckedAt(ftrack))->GetPdgCode(); |
2625 | |
2626 | if (iparent >= 0) { |
2627 | Int_t ip; |
2628 | while(1) { |
2629 | ip=((TParticle*) fPartArray-> |
2630 | UncheckedAt(iparent))->GetFirstMother(); |
2631 | if (ip < 0) { |
2632 | id2 = ((TParticle*) fPartArray-> |
2633 | UncheckedAt(iparent))->GetPdgCode(); |
2634 | break; |
2635 | } else { |
2636 | iparent = ip; |
2637 | id2 = ((TParticle*) fPartArray-> |
2638 | UncheckedAt(iparent))->GetPdgCode(); |
2639 | } // ip<0 |
2640 | } // while |
2641 | }// iparent |
2642 | Int_t id1 = ftrack; |
2643 | Int_t idum = trsign+1; |
2644 | |
2645 | if (id2==443 || id2==553) { |
2646 | nres++; |
2647 | if (id2==443) id2=114; |
2648 | else id2=116; |
2649 | } |
2650 | |
2651 | idg[maxidg] = 30000*id1+10000*idum+id2; |
2652 | pvert1g[maxidg] = (Double_t) Part->Py(); |
2653 | pvert2g[maxidg] = (Double_t) Part->Px(); |
2654 | pvert3g[maxidg] = (Double_t) Part->Pz(); |
2655 | zvertg[maxidg] = (Double_t) Part->Vz(); |
2656 | } //if muon |
2657 | } //if nch |
2658 | maxidg++; |
2659 | } // check limits |
2660 | } // sign+bgr, highest bgr |
2661 | } |
2662 | //pure resonance or mixed cluster with the highest |
2663 | //contribution coming from resonance |
2664 | if (mRaw->fClusterType >= 1 && ftrack1>=0) |
2665 | { |
2666 | if (maxidg<=20000) { |
2667 | if (flag == 0) { |
2668 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2669 | ytrg[maxidg] = (Double_t) mCor->fX[0]; |
2670 | }else { |
2671 | xtrg[maxidg] = (Double_t) mCor->fY[3]; |
2672 | ytrg[maxidg] = (Double_t) mCor->fX[3]; |
2673 | } |
2674 | clsize1[maxidg] = mult1; |
2675 | clsize2[maxidg] = mult2; |
2676 | izch[maxidg] = ich+1; |
2677 | |
2678 | Int_t nent=(Int_t)gAlice->TreeD()->GetEntries(); |
2679 | gAlice->ResetDigits(); |
2680 | if (flag==0) { |
2681 | //gAlice->TreeD()->GetEvent(2); // cathode 2 |
2682 | gAlice->TreeD()->GetEvent(nent-1); // cathode 2 |
2683 | } else { |
2684 | //gAlice->TreeD()->GetEvent(1); // cathode 1 |
2685 | gAlice->TreeD()->GetEvent(nent-2); // cathode 1 |
2686 | } |
2687 | |
2688 | TClonesArray *MUONdigits = MUON->DigitsAddress(ich); |
2689 | Int_t mul=mRaw->fMultiplicity; |
2690 | for (int i=0;i<mul;i++) { |
2691 | Int_t idx=mRaw->fIndexMap[i]; |
2692 | AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx); |
2693 | ihit=dig->fHit-1; |
2694 | if (ihit >= 0) break; |
2695 | |
2696 | } // loop over indexmap |
2697 | //printf("fClusterType, ihit %d %d \n",mRaw->fClusterType,ihit); |
2698 | if (ihit < 0) { |
2699 | xgeant[maxidg] = 0; // only for resonances |
2700 | ygeant[maxidg] = 0; // only for resonances |
2701 | |
2702 | cx[maxidg] = 0; // Px/P of hit |
2703 | cy[maxidg] = 0; // Py/P of hit |
2704 | cz[maxidg] = 0; // Pz/P of hit |
2705 | itypg[maxidg] = 0; // particle number |
2706 | ptotg[maxidg] = 0; // P of hit |
2707 | idg[maxidg] = 0; |
2708 | pvert1g[maxidg] = 0; // Px vertex |
2709 | pvert2g[maxidg] = 0; // Py vertex |
2710 | pvert3g[maxidg] = 0; // Pz vertex |
2711 | zvertg[maxidg] = 0; |
2712 | } else { |
2713 | gAlice->ResetHits(); |
2714 | TH->GetEvent(ftrack1); |
2715 | TClonesArray *MUONhits = MUON->Hits(); |
2716 | AliMUONhit *mHit= (AliMUONhit*)MUONhits-> |
2717 | UncheckedAt(ihit); |
2718 | TClonesArray *fPartArray = gAlice->Particles(); |
2719 | TParticle *Part; |
2720 | Int_t nch=mHit->fChamber-1; |
2721 | //printf("signal ich, nch %d %d \n",ich,nch); |
2722 | if (nch==ich) { |
2723 | Int_t ftrack = mHit->fTrack; |
2724 | Int_t id = ((TParticle*) fPartArray-> |
2725 | UncheckedAt(ftrack))->GetPdgCode(); |
2726 | //printf("id %d \n",id); |
2727 | if (id==kMuonPlus||id==kMuonMinus) { |
2728 | xgeant[maxidg] = (Double_t) mHit->fY; |
2729 | ygeant[maxidg] = (Double_t) mHit->fX; |
2730 | cx[maxidg] = (Double_t) mHit->fCyHit; |
2731 | cy[maxidg] = (Double_t) mHit->fCxHit; |
2732 | cz[maxidg] = (Double_t) mHit->fCzHit; |
2733 | |
2734 | if (id==kMuonPlus) { |
2735 | itypg[maxidg] = 5; |
2736 | } else if (id==kMuonMinus) { |
2737 | itypg[maxidg] = 6; |
2738 | } else itypg[maxidg] = 0; |
2739 | |
2740 | ptotg[maxidg] = (Double_t) mHit->fPTot; |
2741 | Part = (TParticle*) fPartArray-> |
2742 | UncheckedAt(ftrack); |
2743 | Int_t iparent = Part->GetFirstMother(); |
2744 | Int_t id2; |
2745 | id2 = ((TParticle*) fPartArray-> |
2746 | UncheckedAt(ftrack))->GetPdgCode(); |
2747 | |
2748 | if (iparent >= 0) { |
2749 | Int_t ip; |
2750 | while(1) { |
2751 | ip=((TParticle*) fPartArray-> |
2752 | UncheckedAt(iparent))->GetFirstMother(); |
2753 | if (ip < 0) { |
2754 | id2 = ((TParticle*) fPartArray-> |
2755 | UncheckedAt(iparent))->GetPdgCode(); |
2756 | break; |
2757 | } else { |
2758 | iparent = ip; |
2759 | id2 = ((TParticle*) fPartArray-> |
2760 | UncheckedAt(iparent))->GetPdgCode(); |
2761 | } // ip<0 |
2762 | } // while |
2763 | }// iparent |
2764 | Int_t id1 = ftrack; |
2765 | Int_t idum = ftrack1+1; |
2766 | |
2767 | if (id2==443 || id2==553) { |
2768 | nres++; |
2769 | if (id2==443) id2=114; |
2770 | else id2=116; |
2771 | } |
2772 | // printf("id2 %d\n",id2); |
2773 | idg[maxidg] = 30000*id1+10000*idum+id2; |
2774 | pvert1g[maxidg] = (Double_t) Part->Py(); |
2775 | pvert2g[maxidg] = (Double_t) Part->Px(); |
2776 | pvert3g[maxidg] = (Double_t) Part->Pz(); |
2777 | zvertg[maxidg] = (Double_t) Part->Vz(); |
2778 | } //if muon |
2779 | } //if nch |
2780 | } // ihit |
2781 | maxidg++; |
2782 | } // check limits |
2783 | } // if cluster type |
2784 | } // icor loop |
2785 | } // ich loop |
2786 | }// if MUON |
2787 | |
2788 | |
2789 | ievr = nev; |
2790 | cout<<"evenement "<<ievr<<endl; |
2791 | nhittot1 = maxidg ; |
2792 | cout<<"nhittot1="<<nhittot1<<endl; |
2793 | |
2794 | static Int_t nbres=0; |
2795 | static Int_t nbcor=0; |
2796 | if (nres>=19) nbres++; |
2797 | printf("nres ,nncor - %d %d\n",nres,nncor); |
2798 | printf("nbres - %d\n",nbres); |
2799 | if (nncor>=20) nbcor++; |
2800 | printf("nbcor - %d\n",nbcor); |
2801 | printf("nuncor - %d\n",nuncor); |
2802 | printf("nbadcor - %d\n",nbadcor); |
2803 | |
2804 | TC->Reset(); |
2805 | |
2806 | hfile_global->cd(); |
2807 | |
2808 | } |
2809 | |
2810 | void trackf_fit(Int_t &ivertex, Double_t *pest, Double_t *pstep, Double_t &pxzinv, Double_t &tphi, Double_t &talam, Double_t &xvert, Double_t &yvert) |
2811 | { |
2812 | // |
2813 | // Fit a track candidate with the following input parameters: |
2814 | // INPUT : IVERTEX : vertex flag, if IVERTEX=1 (XVERT,YVERT) are free paramaters |
2815 | // if IVERTEX=1 (XVERT,YVERT)=(0.,0.) |
2816 | // PEST(5) : starting value of parameters (minuit) |
2817 | // PSTEP(5) : step size for parameters (minuit) |
2818 | // OUTPUT : PXZINV,TPHI,TALAM,XVERT,YVERT : fitted value of the parameters |
2819 | |
2820 | static Double_t arglist[10]; |
2821 | static Double_t c[5] = {0.4, 0.45, 0.45, 90., 90.}; |
2822 | static Double_t b1, b2, epxz, efi, exs, exvert, eyvert; |
2823 | TString chname; |
2824 | Int_t ierflg = 0; |
2825 | |
2826 | TMinuit *gMinuit = new TMinuit(5); |
2827 | gMinuit->mninit(5,10,7); |
2828 | gMinuit->SetFCN(fcnf); // constant m.f. |
2829 | |
2830 | arglist[0] = -1; |
2831 | |
2832 | gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg); |
2833 | // gMinuit->mnseti('track fitting'); |
2834 | |
2835 | gMinuit->mnparm(0, "invmom", pest[0], pstep[0], -c[0], c[0], ierflg); |
2836 | gMinuit->mnparm(1, "azimuth", pest[1], pstep[1], -c[1], c[1], ierflg); |
2837 | gMinuit->mnparm(2, "deep", pest[2], pstep[2], -c[2], c[2], ierflg); |
2838 | if (ivertex==1) { |
2839 | gMinuit->mnparm(3, "x ", pest[3], pstep[3], -c[3], c[3], ierflg); |
2840 | gMinuit->mnparm(4, "y ", pest[4], pstep[4], -c[4], c[4], ierflg); |
2841 | } |
2842 | |
2843 | gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg); |
2844 | gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg); |
2845 | gMinuit->mnexcm("EXIT" , arglist, 0, ierflg); |
2846 | |
2847 | gMinuit->mnpout(0, chname, pxzinv, epxz, b1, b2, ierflg); |
2848 | gMinuit->mnpout(1, chname, tphi, efi, b1, b2, ierflg); |
2849 | gMinuit->mnpout(2, chname, talam, exs, b1, b2, ierflg); |
2850 | if (ivertex==1) { |
2851 | gMinuit->mnpout(3, chname, xvert, exvert, b1, b2, ierflg); |
2852 | gMinuit->mnpout(4, chname, yvert, eyvert, b1, b2, ierflg); |
2853 | } |
2854 | |
2855 | delete gMinuit; |
2856 | |
2857 | } |
2858 | |
2859 | void fcnf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *pest, Int_t iflag) |
2860 | { |
2861 | // |
2862 | // function called by trackf_fit |
2863 | Int_t futil = 0; |
2864 | fcn(npar,grad,fval,pest,iflag,futil); |
2865 | } |
2866 | |
2867 | void prec_fit(Double_t &pxzinv, Double_t &fis, Double_t &alams, Double_t &xvert, Double_t &yvert, Double_t &pxzinvf, Double_t &fif, Double_t &alf, Double_t &xvertf, Double_t &yvertf, Double_t &epxzinv, Double_t &efi, Double_t &exs, Double_t &exvert, Double_t &eyvert) |
2868 | { |
2869 | // |
2870 | // minuit fits for tracking finding |
2871 | |
2872 | static Double_t arglist[10]; |
2873 | static Double_t c1[5] = {0.001, 0.001, 0.001, 1., 1.}; |
2874 | static Double_t c2[5] = {0.5, 0.5, 0.5, 120., 120.}; |
2875 | static Double_t emat[9]; |
2876 | static Double_t b1, b2; |
2877 | Double_t fmin, fedm, errdef; |
2878 | Int_t npari, nparx, istat; |
2879 | |
2880 | TString chname; |
2881 | Int_t ierflg = 0; |
2882 | |
2883 | TMinuit *gMinuit = new TMinuit(5); |
2884 | gMinuit->mninit(5,10,7); |
2885 | gMinuit->SetFCN(fcnfitf); |
2886 | |
2887 | arglist[0] = -1.; |
2888 | gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg); |
2889 | |
2890 | // gMinuit->mnseti('track fitting'); |
2891 | |
2892 | gMinuit->mnparm(0,"invmom", pxzinv, c1[0], -c2[0], c2[0], ierflg); // 0.003, 0.5 |
2893 | gMinuit->mnparm(1,"azimuth ", fis, c1[1], -c2[1], c2[1], ierflg); |
2894 | gMinuit->mnparm(2,"deep ", alams, c1[2], -c2[2], c2[2], ierflg); |
2895 | gMinuit->mnparm(3,"xvert", xvert, c1[3], -c2[3], c2[3], ierflg); |
2896 | gMinuit->mnparm(4,"yvert", yvert, c1[4], -c2[4], c2[4], ierflg); |
2897 | |
2898 | gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg); |
2899 | arglist[0] = 2.; |
2900 | gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg); |
2901 | gMinuit->mnexcm("EXIT", arglist, 0, ierflg); |
2902 | |
2903 | gMinuit->mnpout(0, chname, pxzinvf, epxzinv, b1, b2, ierflg); |
2904 | gMinuit->mnpout(1, chname, fif, efi, b1, b2, ierflg); |
2905 | gMinuit->mnpout(2, chname, alf, exs, b1, b2, ierflg); |
2906 | gMinuit->mnpout(3, chname, xvertf, exvert, b1, b2, ierflg); |
2907 | gMinuit->mnpout(4, chname, yvertf, eyvert, b1, b2, ierflg); |
2908 | |
2909 | gMinuit->mnemat(emat, 3); |
2910 | gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat); |
2911 | |
2912 | delete gMinuit; |
2913 | } |
2914 | |
2915 | void fcnfitf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *xval, Int_t iflag) |
2916 | { |
2917 | // |
2918 | // function called by prec_fit |
2919 | Int_t futil = 0; |
2920 | fcnfit(npar,grad,fval,xval,iflag,futil); |
2921 | } |
2922 | |
2923 | ///////////////////// fin modifs perso ////////////////////// |
2924 | |
fe4da5cc |
2925 | ClassImp(AliMUONcluster) |
2926 | |
2927 | //___________________________________________ |
2928 | AliMUONcluster::AliMUONcluster(Int_t *clhits) |
2929 | { |
2930 | fHitNumber=clhits[0]; |
2931 | fCathode=clhits[1]; |
2932 | fQ=clhits[2]; |
2933 | fPadX=clhits[3]; |
2934 | fPadY=clhits[4]; |
2935 | fQpad=clhits[5]; |
2936 | fRSec=clhits[6]; |
2937 | } |
2938 | ClassImp(AliMUONdigit) |
2939 | //_____________________________________________________________________________ |
2940 | AliMUONdigit::AliMUONdigit(Int_t *digits) |
2941 | { |
2942 | // |
2943 | // Creates a MUON digit object to be updated |
2944 | // |
a897a37a |
2945 | fPadX = digits[0]; |
2946 | fPadY = digits[1]; |
2947 | fSignal = digits[2]; |
2948 | fPhysics = digits[3]; |
2949 | fHit = digits[4]; |
fe4da5cc |
2950 | |
2951 | } |
2952 | //_____________________________________________________________________________ |
2953 | AliMUONdigit::AliMUONdigit(Int_t *tracks, Int_t *charges, Int_t *digits) |
2954 | { |
2955 | // |
2956 | // Creates a MUON digit object |
2957 | // |
2958 | fPadX = digits[0]; |
2959 | fPadY = digits[1]; |
2960 | fSignal = digits[2]; |
a897a37a |
2961 | fPhysics = digits[3]; |
2962 | fHit = digits[4]; |
fe4da5cc |
2963 | for(Int_t i=0; i<10; i++) { |
2964 | fTcharges[i] = charges[i]; |
2965 | fTracks[i] = tracks[i]; |
2966 | } |
2967 | } |
2968 | |
a897a37a |
2969 | AliMUONdigit::~AliMUONdigit() |
2970 | { |
2971 | |
2972 | } |
2973 | |
fe4da5cc |
2974 | ClassImp(AliMUONlist) |
2975 | |
2976 | //____________________________________________________________________________ |
a897a37a |
2977 | AliMUONlist::AliMUONlist(Int_t ich, Int_t *digits): |
fe4da5cc |
2978 | AliMUONdigit(digits) |
2979 | { |
2980 | // |
2981 | // Creates a MUON digit list object |
2982 | // |
2983 | |
a897a37a |
2984 | fChamber = ich; |
fe4da5cc |
2985 | fTrackList = new TObjArray; |
2986 | |
2987 | } |
fe4da5cc |
2988 | |
2989 | ClassImp(AliMUONhit) |
2990 | |
2991 | //___________________________________________ |
2992 | AliMUONhit::AliMUONhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits): |
2993 | AliHit(shunt, track) |
2994 | { |
2995 | fChamber=vol[0]; |
a897a37a |
2996 | fParticle=hits[0]; |
fe4da5cc |
2997 | fX=hits[1]; |
2998 | fY=hits[2]; |
2999 | fZ=hits[3]; |
3000 | fTheta=hits[4]; |
3001 | fPhi=hits[5]; |
3002 | fTlength=hits[6]; |
3003 | fEloss=hits[7]; |
3004 | fPHfirst=(Int_t) hits[8]; |
3005 | fPHlast=(Int_t) hits[9]; |
fe4da5cc |
3006 | |
a897a37a |
3007 | // modifs perso |
3008 | fPTot=hits[10]; |
3009 | fCxHit=hits[11]; |
3010 | fCyHit=hits[12]; |
3011 | fCzHit=hits[13]; |
3012 | } |
3013 | ClassImp(AliMUONcorrelation) |
3014 | //___________________________________________ |
3015 | //_____________________________________________________________________________ |
3016 | AliMUONcorrelation::AliMUONcorrelation(Int_t *idx, Float_t *x, Float_t *y) |
3017 | { |
3018 | // |
3019 | // Creates a MUON correlation object |
3020 | // |
3021 | for(Int_t i=0; i<4; i++) { |
3022 | fCorrelIndex[i] = idx[i]; |
3023 | fX[i] = x[i]; |
3024 | fY[i] = y[i]; |
3025 | } |
3026 | } |
3027 | ClassImp(AliMUONRawCluster) |
3028 | Int_t AliMUONRawCluster::Compare(TObject *obj) |
fe4da5cc |
3029 | { |
a897a37a |
3030 | /* |
3031 | AliMUONRawCluster *raw=(AliMUONRawCluster *)obj; |
3032 | Float_t r=GetRadius(); |
3033 | Float_t ro=raw->GetRadius(); |
3034 | if (r>ro) return 1; |
3035 | else if (r<ro) return -1; |
3036 | else return 0; |
3037 | */ |
3038 | AliMUONRawCluster *raw=(AliMUONRawCluster *)obj; |
3039 | Float_t y=fY; |
3040 | Float_t yo=raw->fY; |
3041 | if (y>yo) return 1; |
3042 | else if (y<yo) return -1; |
3043 | else return 0; |
3044 | |
fe4da5cc |
3045 | } |
3046 | |
a897a37a |
3047 | Int_t AliMUONRawCluster:: |
3048 | BinarySearch(Float_t y, TArrayF coord, Int_t from, Int_t upto) |
fe4da5cc |
3049 | { |
a897a37a |
3050 | // Find object using a binary search. Array must first have been sorted. |
3051 | // Search can be limited by setting upto to desired index. |
3052 | |
3053 | Int_t low=from, high=upto-1, half; |
3054 | while(high-low>1) { |
3055 | half=(high+low)/2; |
3056 | if(y>coord[half]) low=half; |
3057 | else high=half; |
3058 | } |
3059 | return low; |
fe4da5cc |
3060 | } |
3061 | |
a897a37a |
3062 | void AliMUONRawCluster::SortMin(Int_t *idx,Float_t *xdarray,Float_t *xarray,Float_t *yarray,Float_t *qarray, Int_t ntr) |
fe4da5cc |
3063 | { |
a897a37a |
3064 | // |
3065 | // Get the 3 closest points(cog) one can find on the second cathode |
3066 | // starting from a given cog on first cathode |
3067 | // |
3068 | |
3069 | // |
3070 | // Loop over deltax, only 3 times |
3071 | // |
3072 | |
3073 | Float_t xmin; |
3074 | Int_t jmin; |
3075 | Int_t id[3] = {-2,-2,-2}; |
3076 | Float_t jx[3] = {0.,0.,0.}; |
3077 | Float_t jy[3] = {0.,0.,0.}; |
3078 | Float_t jq[3] = {0.,0.,0.}; |
3079 | Int_t jid[3] = {-2,-2,-2}; |
3080 | Int_t i,j,imax; |
3081 | |
3082 | if (ntr<3) imax=ntr; |
3083 | else imax=3; |
3084 | for(i=0;i<imax;i++){ |
3085 | xmin=1001.; |
3086 | jmin=0; |
3087 | |
3088 | for(j=0;j<ntr;j++){ |
3089 | if ((i == 1 && j == id[i-1]) |
3090 | ||(i == 2 && (j == id[i-1] || j == id[i-2]))) continue; |
3091 | if (TMath::Abs(xdarray[j]) < xmin) { |
3092 | xmin = TMath::Abs(xdarray[j]); |
3093 | jmin=j; |
3094 | } |
3095 | } // j |
3096 | if (xmin != 1001.) { |
3097 | id[i]=jmin; |
3098 | jx[i]=xarray[jmin]; |
3099 | jy[i]=yarray[jmin]; |
3100 | jq[i]=qarray[jmin]; |
3101 | jid[i]=idx[jmin]; |
3102 | } |
3103 | |
3104 | } // i |
3105 | |
3106 | for (i=0;i<3;i++){ |
3107 | if (jid[i] == -2) { |
3108 | xarray[i]=1001.; |
3109 | yarray[i]=1001.; |
3110 | qarray[i]=1001.; |
3111 | idx[i]=-1; |
3112 | } else { |
3113 | xarray[i]=jx[i]; |
3114 | yarray[i]=jy[i]; |
3115 | qarray[i]=jq[i]; |
3116 | idx[i]=jid[i]; |
3117 | } |
fe4da5cc |
3118 | } |
a897a37a |
3119 | |
fe4da5cc |
3120 | } |
3121 | |
3122 | |
a897a37a |
3123 | Int_t AliMUONRawCluster::PhysicsContribution() |
fe4da5cc |
3124 | { |
a897a37a |
3125 | Int_t iPhys=0; |
3126 | Int_t iBg=0; |
3127 | Int_t iMixed=0; |
3128 | for (Int_t i=0; i<fMultiplicity; i++) { |
3129 | if (fPhysicsMap[i]==2) iPhys++; |
3130 | if (fPhysicsMap[i]==1) iMixed++; |
3131 | if (fPhysicsMap[i]==0) iBg++; |
3132 | } |
3133 | if (iMixed==0 && iBg==0) { |
3134 | return 2; |
3135 | } else if ((iPhys != 0 && iBg !=0) || iMixed != 0) { |
3136 | return 1; |
3137 | } else { |
3138 | return 0; |
3139 | } |
fe4da5cc |
3140 | } |
3141 | |
a897a37a |
3142 | |
3143 | ClassImp(AliMUONreccluster) |
3144 | ClassImp(AliMUONsegmentation) |
3145 | ClassImp(AliMUONresponse) |
3146 | |
3147 | |
3148 | |
3149 | |
3150 | |
fe4da5cc |
3151 | |
fe4da5cc |
3152 | |
fe4da5cc |
3153 | |
fe4da5cc |
3154 | |
3155 | |
3156 | |
3157 | |
3158 | |
3159 | |
3160 | |
3161 | |
3162 | |
3163 | |
3164 | |