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