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