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