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