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