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