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