1 ////////////////////////////////////////////////
2 // Manager and hits classes for set:MUON //
3 ////////////////////////////////////////////////
7 #include <TRotMatrix.h>
13 #include <TObjArray.h>
15 #include <TParticle.h>
21 #include <TDirectory.h>
22 #include <TObjectTable.h>
27 #include "AliMUONClusterFinder.h"
31 #include "AliCallf77.h"
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_
46 # define hist_create hist_create_
47 # define hist_closed hist_closed_
50 # define trackf_fit trackf_fit_
51 # define prec_fit prec_fit_
52 # define fcnfit fcnfit_
53 # define reco_term reco_term_
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
67 # define hist_create HIST_CREATE
68 # define hist_closed HIST_CLOSED
71 # define trackf_fit TRACKF_FIT
72 # define prec_fit PREC_FIT
73 # define fcnfit FCNFIT
74 # define reco_term RECO_TERM
79 void type_of_call reco_init(Double_t &, Double_t &, Double_t &);
80 void type_of_call reco_term();
81 void type_of_call cutpxz(Double_t &);
82 void type_of_call sigmacut(Double_t &);
83 void type_of_call xpreci(Double_t &);
84 void type_of_call ypreci(Double_t &);
85 void type_of_call reconstmuon(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
86 void 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 *);
87 void 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 *);
88 void type_of_call chfill(Int_t &, Float_t &, Float_t &, Float_t &);
89 void type_of_call chfill2(Int_t &, Float_t &, Float_t &, Float_t &);
90 void type_of_call chf1(Int_t &, Float_t &, Float_t &);
91 void 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 *);
92 void type_of_call hist_create();
93 void type_of_call hist_closed();
94 void type_of_call fcnf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
95 void type_of_call fcn(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &);
96 void type_of_call trackf_fit(Int_t &, Double_t *, Double_t *, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
97 void 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 &);
98 void type_of_call fcnfitf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
99 void type_of_call fcnfit(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &);
100 Float_t type_of_call rndm() {return gRandom->Rndm();}
103 // Static variables for the pad-hit iterator routines
104 static Int_t sMaxIterPad=0;
105 static Int_t sCurIterPad=0;
108 static TClonesArray *fHits2; //Listof hits for one track only
109 static TClonesArray *fClusters2; //List of clusters for one track only
110 static TClonesArray *fParticles2; //List of particles in the Kine tree
112 //___________________________________________
140 //___________________________________________
141 AliMUON::AliMUON(const char *name, const char *title)
142 : AliDetector(name,title)
146 <img src="gif/alimuon.gif">
150 fHits = new TClonesArray("AliMUONhit",1000);
151 fClusters = new TClonesArray("AliMUONcluster",10000);
155 fNdch = new Int_t[10];
157 fDchambers = new TObjArray(10);
161 for (i=0; i<10 ;i++) {
162 (*fDchambers)[i] = new TClonesArray("AliMUONdigit",10000);
166 fNrawch = new Int_t[10];
168 fRawClusters = new TObjArray(10);
170 for (i=0; i<10 ;i++) {
171 (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000);
175 fNcorch = new Int_t[10];
176 fCathCorrel = new TObjArray(10);
177 for (i=0; i<10 ;i++) {
178 (*fCathCorrel)[i] = new TClonesArray("AliMUONcorrelation",1000);
185 // Transport angular cut
200 SetMarkerColor(kRed);
203 //___________________________________________
207 printf("Calling AliMUON destructor !!!\n");
216 delete (*fDchambers)[i];
222 delete (*fRawClusters)[i];
228 delete (*fCathCorrel)[i];
234 //___________________________________________
235 void AliMUON::AddHit(Int_t track, Int_t *vol, Float_t *hits)
237 TClonesArray &lhits = *fHits;
238 new(lhits[fNhits++]) AliMUONhit(fIshunt,track,vol,hits);
240 //___________________________________________
241 void AliMUON::AddCluster(Int_t *clhits)
243 TClonesArray &lclusters = *fClusters;
244 new(lclusters[fNclusters++]) AliMUONcluster(clhits);
246 //_____________________________________________________________________________
247 void AliMUON::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
250 // Add a MUON digit to the list
253 TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
254 new(ldigits[fNdch[id]++]) AliMUONdigit(tracks,charges,digits);
257 //_____________________________________________________________________________
258 void AliMUON::AddRawCluster(Int_t id, const AliMUONRawCluster& c)
261 // Add a MUON digit to the list
264 TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
265 new(lrawcl[fNrawch[id]++]) AliMUONRawCluster(c);
267 //_____________________________________________________________________________
268 void AliMUON::AddCathCorrel(Int_t id, Int_t *idx, Float_t *x, Float_t *y)
271 // Add a MUON digit to the list
274 TClonesArray &lcorrel = *((TClonesArray*)(*fCathCorrel)[id]);
275 new(lcorrel[fNcorch[id]++]) AliMUONcorrelation(idx,x,y);
279 //_____________________________________________________________________________
280 void AliMUON::AddRecCluster(Int_t iCh, Int_t iCat, AliMUONRecCluster* Cluster)
283 // Add a MUON reconstructed cluster to the list
286 TObjArray* ClusterList = RecClusters(iCh,iCat);
287 ClusterList->Add(Cluster);
291 //___________________________________________
292 void AliMUON::BuildGeometry()
294 TNode *Node, *NodeF, *Top;
295 const int kColorMUON = kBlue;
297 Top=gAlice->GetGeometry()->GetNode("alice");
300 // z-Positions of Chambers
301 const Float_t cz[5]={511., 686., 971., 1245., 1445.};
304 const Float_t dmi[5]={ 35., 47., 67., 86., 100.};
307 const Float_t dma[5]={183., 245., 346., 520., 520.};
309 TRotMatrix* rot000 = new TRotMatrix("Rot000"," ", 90, 0, 90, 90, 0, 0);
310 TRotMatrix* rot090 = new TRotMatrix("Rot090"," ", 90, 90, 90,180, 0, 0);
311 TRotMatrix* rot180 = new TRotMatrix("Rot180"," ", 90,180, 90,270, 0, 0);
312 TRotMatrix* rot270 = new TRotMatrix("Rot270"," ", 90,270, 90, 0, 0, 0);
315 float rmin, rmax, dx, dy, dz, dr, zpos;
317 char NameChamber[9], NameSense[9], NameFrame[9], NameNode[7];
318 for (Int_t i=0; i<5; i++) {
319 for (Int_t j=0; j<2; j++) {
328 sprintf(NameChamber,"C_MUON%d",id);
329 sprintf(NameSense,"S_MUON%d",id);
330 sprintf(NameFrame,"F_MUON%d",id);
333 new TTUBE(NameChamber,"Mother","void",rmin,rmax,0.25,1.);
336 new TTUBE(NameSense,"Sens. region","void",rmin,rmax,0.25, 1.);
340 TBRIK* FMUON = new TBRIK(NameFrame,"Frame","void",dx,dy,dz);
342 sprintf(NameNode,"MUON%d",100+id);
343 Node = new TNode(NameNode,"ChamberNode",NameChamber,0,0,zpos,"");
344 Node->SetLineColor(kColorMUON);
347 sprintf(NameNode,"MUON%d",200+id);
348 Node = new TNode(NameNode,"Sens. Region Node",NameSense,0,0,0,"");
349 Node->SetLineColor(kColorMUON);
352 sprintf(NameNode,"MUON%d",300+id);
353 NodeF = new TNode(NameNode,"Frame0",FMUON,dr, 0, 0,rot000,"");
354 NodeF->SetLineColor(kColorMUON);
356 sprintf(NameNode,"MUON%d",400+id);
357 NodeF = new TNode(NameNode,"Frame1",FMUON,0 ,dr,0,rot090,"");
358 NodeF->SetLineColor(kColorMUON);
360 sprintf(NameNode,"MUON%d",500+id);
361 NodeF = new TNode(NameNode,"Frame2",FMUON,-dr,0,0,rot180,"");
362 NodeF->SetLineColor(kColorMUON);
364 sprintf(NameNode,"MUON%d",600+id);
365 NodeF = new TNode(NameNode,"Frame3",FMUON,0,-dr,0,rot270,"");
366 NodeF->SetLineColor(kColorMUON);
372 //___________________________________________
373 Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
378 //___________________________________________
379 void AliMUON::MakeBranch(Option_t* option)
381 // Create Tree branches for the MUON.
383 const Int_t buffersize = 4000;
385 sprintf(branchname,"%sCluster",GetName());
387 AliDetector::MakeBranch(option);
389 if (fClusters && gAlice->TreeH()) {
390 gAlice->TreeH()->Branch(branchname,&fClusters, buffersize);
391 printf("Making Branch %s for clusters\n",branchname);
394 // one branch for digits per chamber
397 for (i=0; i<10 ;i++) {
398 sprintf(branchname,"%sDigits%d",GetName(),i+1);
400 if (fDchambers && gAlice->TreeD()) {
401 gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), buffersize);
402 printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
406 printf("Make Branch - TreeR address %p\n",gAlice->TreeR());
408 // one branch for raw clusters per chamber
409 for (i=0; i<10 ;i++) {
410 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
412 if (fRawClusters && gAlice->TreeR()) {
413 gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), buffersize);
414 printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
418 // Emmanuel's stuff - one branch for rec clusters
420 for (i=0; i<20; i++) {
421 sprintf(branchname,"%sRecClus%d",GetName(),i+1);
422 if (fRecClusters && gAlice->TreeD()) {
424 ->Branch(branchname,"TObjArray",
425 &((*fRecClusters)[i]), buffersize,0);
426 printf("Making Branch %s for clusters in chamber %d\n",
434 //___________________________________________
435 void AliMUON::SetTreeAddress()
437 // Set branch address for the Hits and Digits Tree.
439 AliDetector::SetTreeAddress();
442 TTree *treeH = gAlice->TreeH();
443 TTree *treeD = gAlice->TreeD();
444 TTree *treeR = gAlice->TreeR();
448 branch = treeH->GetBranch("MUONCluster");
449 if (branch) branch->SetAddress(&fClusters);
454 for (int i=0; i<10; i++) {
455 sprintf(branchname,"%sDigits%d",GetName(),i+1);
457 branch = treeD->GetBranch(branchname);
458 if (branch) branch->SetAddress(&((*fDchambers)[i]));
463 // printf("SetTreeAddress --- treeR address %p \n",treeR);
466 for (int i=0; i<10; i++) {
467 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
469 branch = treeR->GetBranch(branchname);
470 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
478 for (int i=0; i<20; i++) {
479 sprintf(branchname,"%sRecClus%d",GetName(),i+1);
481 branch = treeR->GetBranch(branchname);
482 if (branch) branch->SetAddress(&((*fRecClusters)[i]));
488 //___________________________________________
489 void AliMUON::ResetHits()
491 // Reset number of clusters and the cluster array for this detector
492 AliDetector::ResetHits();
494 if (fClusters) fClusters->Clear();
497 //____________________________________________
498 void AliMUON::ResetDigits()
501 // Reset number of digits and the digits array for this detector
503 for ( int i=0;i<10;i++ ) {
504 if ((*fDchambers)[i]) ((TClonesArray*)(*fDchambers)[i])->Clear();
505 if (fNdch) fNdch[i]=0;
508 //____________________________________________
509 void AliMUON::ResetRawClusters()
512 // Reset number of raw clusters and the raw clust array for this detector
514 for ( int i=0;i<10;i++ ) {
515 if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
516 if (fNrawch) fNrawch[i]=0;
519 //____________________________________________
520 void AliMUON::ResetCorrelation()
523 // Reset number of correl clusters and the correl clust array for
526 for ( int i=0;i<10;i++ ) {
527 if ((*fCathCorrel)[i]) ((TClonesArray*)(*fCathCorrel)[i])->Clear();
528 if (fNcorch) fNcorch[i]=0;
532 //____________________________________________
533 void AliMUON::ResetRecClusters()
536 // Reset the rec clusters
539 for ( int i=0;i<20;i++ ) {
540 if ((*fRecClusters)[i]) (*fRecClusters)[i]->Clear();
544 //___________________________________________
546 void AliMUON::SetPADSIZ(Int_t id, Int_t isec, Float_t p1, Float_t p2)
549 ((AliMUONchamber*) (*fChambers)[i]) ->SetPADSIZ(isec,p1,p2);
550 ((AliMUONchamber*) (*fChambers)[i+1])->SetPADSIZ(isec,p1,p2);
553 //___________________________________________
554 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
557 ((AliMUONchamber*) (*fChambers)[i])->SetChargeSlope(p1);
558 ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSlope(p1);
561 //___________________________________________
562 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
565 ((AliMUONchamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
566 ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSpread(p1,p2);
569 //___________________________________________
570 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
573 ((AliMUONchamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
574 ((AliMUONchamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
577 //___________________________________________
578 void AliMUON::SetMaxAdc(Int_t id, Float_t p1)
581 ((AliMUONchamber*) (*fChambers)[i])->SetMaxAdc(p1);
582 ((AliMUONchamber*) (*fChambers)[i+1])->SetMaxAdc(p1);
585 //___________________________________________
586 void AliMUON::SetMaxStepGas(Float_t p1)
591 //___________________________________________
592 void AliMUON::SetMaxStepAlu(Float_t p1)
597 //___________________________________________
598 void AliMUON::SetMaxDestepGas(Float_t p1)
603 //___________________________________________
604 void AliMUON::SetMaxDestepAlu(Float_t p1)
608 //___________________________________________
609 void AliMUON::SetMuonAcc(Bool_t acc, Float_t angmin, Float_t angmax)
615 //___________________________________________
616 void AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliMUONsegmentation *segmentation)
618 ((AliMUONchamber*) (*fChambers)[id])->SegmentationModel(isec, segmentation);
621 //___________________________________________
622 void AliMUON::SetResponseModel(Int_t id, AliMUONresponse *response)
624 ((AliMUONchamber*) (*fChambers)[id])->ResponseModel(response);
627 void AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinder *reconst)
629 ((AliMUONchamber*) (*fChambers)[id])->ReconstructionModel(reconst);
632 void AliMUON::SetNsec(Int_t id, Int_t nsec)
634 ((AliMUONchamber*) (*fChambers)[id])->SetNsec(nsec);
638 //___________________________________________
640 void AliMUON::StepManager()
642 printf("Dummy version of muon step -- it should never happen!!\n");
643 const Float_t kRaddeg = 180/TMath::Pi();
644 AliMC* pMC = AliMC::GetMC();
647 Float_t pt, th0, th2;
650 if((nsec=pMC->NSecondaries())>0) {
651 pMC->ProdProcess(proc);
652 if((pMC->TrackPid()==443 || pMC->TrackPid()==553) && !strcmp(proc,"DCAY")) {
654 // Check angular acceptance
655 //* --- and have muons from resonance decays in the wanted window ---
657 printf(" AliMUON::StepManager: Strange resonance Decay into %d particles\n",nsec);
660 pMC->GetSecondary(0,ipart,x,p);
661 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
662 th0 = TMath::ATan2(pt,p[2])*kRaddeg;
663 pMC->GetSecondary(1,ipart,x,p);
664 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
665 th2 = TMath::ATan2(pt,p[2])*kRaddeg;
666 if(!(fAccMin < th0 && th0 < fAccMax) ||
667 !(fAccMin < th2 && th2 < fAccMax))
675 void AliMUON::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol)
678 // Calls the charge disintegration method of the current chamber and adds
679 // the simulated cluster to the root treee
682 Float_t newclust[6][500];
687 // Integrated pulse height on chamber
693 ((AliMUONchamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust);
694 // printf("\n Add new clusters %d %f \n", nnew, eloss*1.e9);
699 for (Int_t i=0; i<nnew; i++) {
700 if (Int_t(newclust[3][i]) > 0) {
703 clhits[1] = Int_t(newclust[5][i]);
705 clhits[2] = Int_t(newclust[0][i]);
707 clhits[3] = Int_t(newclust[1][i]);
709 clhits[4] = Int_t(newclust[2][i]);
711 clhits[5] = Int_t(newclust[3][i]);
712 // Pad: chamber sector
713 clhits[6] = Int_t(newclust[4][i]);
718 // printf("\n %d new clusters added", ic);
721 void AliMUON::Digitise(Int_t nev,Int_t bgr_ev,Option_t *option,Option_t *opt,Text_t *filename)
723 // keep galice.root for signal and name differently the file for
724 // background when add! otherwise the track info for signal will be lost !
726 static Bool_t first=kTRUE;
727 // static TTree *TrH1;
729 char *Add = strstr(option,"Add");
730 //char *listoftracks = strstr(opt,"listoftracks");
732 AliMUONchamber* iChamber;
733 AliMUONsegmentation* segmentation;
738 TObjArray *list=new TObjArray;
739 static TClonesArray *p_adr=0;
740 if(!p_adr) p_adr=new TClonesArray("TVector",1000);
743 AliMUON *MUON = (AliMUON *) gAlice->GetModule("MUON");
744 AliMUONHitMap * HitMap[10];
745 for (Int_t i=0; i<10; i++) {HitMap[i]=0;}
749 cout<<"filename"<<fFileName<<endl;
750 File=new TFile(fFileName);
751 cout<<"I have opened "<<fFileName<<" file "<<endl;
752 fHits2 = new TClonesArray("AliMUONhit",1000 );
753 fClusters2 = new TClonesArray("AliMUONcluster",10000);
758 // Get Hits Tree header from file
759 if(fHits2) fHits2->Clear();
760 if(fClusters2) fClusters2->Clear();
761 if(TrH1) delete TrH1;
765 sprintf(treeName,"TreeH%d",bgr_ev);
766 TrH1 = (TTree*)gDirectory->Get(treeName);
767 //printf("TrH1 %p of treename %s for event %d \n",TrH1,treeName,bgr_ev);
770 printf("ERROR: cannot find Hits Tree for event:%d\n",bgr_ev);
772 // Set branch addresses
775 sprintf(branchname,"%s",GetName());
776 if (TrH1 && fHits2) {
777 branch = TrH1->GetBranch(branchname);
778 if (branch) branch->SetAddress(&fHits2);
780 if (TrH1 && fClusters2) {
781 branch = TrH1->GetBranch("MUONCluster");
782 if (branch) branch->SetAddress(&fClusters2);
785 //Int_t ntracks1 =(Int_t)TrH1->GetEntries();
786 //printf("background - ntracks1 - %d\n",ntracks1);
789 // loop over cathodes
793 for (int icat=0; icat<2; icat++) {
795 for (Int_t i=0; i<10; i++) {
805 for (Int_t i =0; i<10; i++) {
806 iChamber=(AliMUONchamber*) (*fChambers)[i];
807 if (iChamber->Nsec()==1 && icat==1) {
810 segmentation=iChamber->GetSegmentationModel(icat+1);
812 HitMap[i] = new AliMUONHitMapA1(segmentation, list);
814 //printf("Start loop over tracks \n");
819 TTree *TH = gAlice->TreeH();
820 Int_t ntracks =(Int_t) TH->GetEntries();
821 //printf("signal - ntracks %d\n",ntracks);
822 Int_t nmuon[10]={0,0,0,0,0,0,0,0,0,0};
826 for (Int_t track=0; track<ntracks; track++) {
832 for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
834 mHit=(AliMUONhit*)MUON->NextHit())
836 Int_t nch = mHit->fChamber-1; // chamber number
837 if (nch >9) continue;
838 iChamber = &(MUON->Chamber(nch));
839 Int_t rmin = (Int_t)iChamber->RInner();
840 Int_t rmax = (Int_t)iChamber->ROuter();
844 if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) {
845 xhit[nch][nmuon[nch]]=mHit->fX;
846 yhit[nch][nmuon[nch]]=mHit->fY;
848 if (nmuon[nch] >2) printf("nmuon %d\n",nmuon[nch]);
857 // Loop over pad hits
858 for (AliMUONcluster* mPad=
859 (AliMUONcluster*)MUON->FirstPad(mHit,fClusters);
861 mPad=(AliMUONcluster*)MUON->NextPad(fClusters))
863 Int_t cathode = mPad->fCathode; // cathode number
864 Int_t ipx = mPad->fPadX; // pad number on X
865 Int_t ipy = mPad->fPadY; // pad number on Y
866 Int_t iqpad = Int_t(mPad->fQpad*kScale);// charge per pad
867 // Int_t iqpad = mPad->fQpad; // charge per pad
871 if (cathode != (icat+1)) continue;
872 // fill the info array
874 segmentation=iChamber->GetSegmentationModel(cathode);
875 segmentation->GetPadCxy(ipx,ipy,thex,they);
876 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
877 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
879 new((*p_adr)[countadr++]) TVector(2);
880 TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]);
881 trinfo(0)=(Float_t)track;
882 trinfo(1)=(Float_t)iqpad;
888 if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) {
889 digits[4]=mPad->fHitNumber;
893 // build the list of fired pads and update the info
894 if (!HitMap[nch]->TestHit(ipx, ipy)) {
896 list->AddAtAndExpand(
897 new AliMUONlist(nch,digits),counter);
899 HitMap[nch]->SetHit(ipx, ipy, counter);
901 pdigit=(AliMUONlist*)list->At(list->GetLast());
903 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
904 trlist->Add(&trinfo);
906 pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy);
908 (*pdigit).fSignal+=iqpad;
909 (*pdigit).fPhysics+=iqpad;
910 // update list of tracks
911 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
912 Int_t last_entry=trlist->GetLast();
913 TVector *ptrk_p=(TVector*)trlist->At(last_entry);
914 TVector &ptrk=*ptrk_p;
915 Int_t last_track=Int_t(ptrk(0));
916 Int_t last_charge=Int_t(ptrk(1));
917 if (last_track==track) {
919 trlist->RemoveAt(last_entry);
920 trinfo(0)=last_track;
921 trinfo(1)=last_charge;
922 trlist->AddAt(&trinfo,last_entry);
924 trlist->Add(&trinfo);
926 // check the track list
927 Int_t nptracks=trlist->GetEntriesFast();
929 for (Int_t tr=0;tr<nptracks;tr++) {
930 TVector *pptrk_p=(TVector*)trlist->At(tr);
931 TVector &pptrk=*pptrk_p;
932 trk[tr]=Int_t(pptrk(0));
933 chtrk[tr]=Int_t(pptrk(1));
937 } //end loop over clusters
941 //Int_t nentr1=list->GetEntriesFast();
942 //printf(" \n counter, nentr1 %d %d\n",counter,nentr1);
944 // open the file with background
947 ntracks =(Int_t)TrH1->GetEntries();
948 //printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
949 //printf("background - Start loop over tracks \n");
953 for (Int_t track=0; track<ntracks; track++) {
955 if (fHits2) fHits2->Clear();
956 if (fClusters2) fClusters2->Clear();
958 TrH1->GetEvent(track);
962 for(int i=0;i<fHits2->GetEntriesFast();++i)
964 mHit=(AliMUONhit*) (*fHits2)[i];
965 Int_t nch = mHit->fChamber-1; // chamber number
966 if (nch >9) continue;
967 iChamber = &(MUON->Chamber(nch));
968 Int_t rmin = (Int_t)iChamber->RInner();
969 Int_t rmax = (Int_t)iChamber->ROuter();
970 Float_t xbgr=mHit->fX;
971 Float_t ybgr=mHit->fY;
974 for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) {
975 Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon])
976 +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]);
977 if (dist<100) cond=kTRUE;
982 // Loop over pad hits
983 //for(int j=0;j<fClusters2->GetEntriesFast();++j)
984 for (AliMUONcluster* mPad=
985 (AliMUONcluster*)MUON->FirstPad(mHit,fClusters2);
987 mPad=(AliMUONcluster*)MUON->NextPad(fClusters2))
989 // mPad = (AliMUONcluster*) (*fClusters2)[j];
990 Int_t cathode = mPad->fCathode; // cathode number
991 Int_t ipx = mPad->fPadX; // pad number on X
992 Int_t ipy = mPad->fPadY; // pad number on Y
993 Int_t iqpad = Int_t(mPad->fQpad*kScale);// charge per pad
994 // Int_t iqpad = mPad->fQpad; // charge per pad
996 if (cathode != (icat+1)) continue;
997 // if (!HitMap[nch]->CheckBoundary()) continue;
998 // fill the info array
1000 segmentation=iChamber->GetSegmentationModel(cathode);
1001 segmentation->GetPadCxy(ipx,ipy,thex,they);
1002 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
1003 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
1005 new((*p_adr)[countadr++]) TVector(2);
1006 TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]);
1007 trinfo(0)=-1; // tag background
1016 AliMUONlist* pdigit;
1017 // build the list of fired pads and update the info
1018 if (!HitMap[nch]->TestHit(ipx, ipy)) {
1019 list->AddAtAndExpand(new AliMUONlist(nch,digits),counter);
1021 HitMap[nch]->SetHit(ipx, ipy, counter);
1024 pdigit=(AliMUONlist*)list->At(list->GetLast());
1026 TObjArray *trlist=(TObjArray*)pdigit->
1028 trlist->Add(&trinfo);
1030 pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy);
1032 (*pdigit).fSignal+=iqpad;
1034 // update list of tracks
1035 TObjArray* trlist=(TObjArray*)pdigit->
1037 Int_t last_entry=trlist->GetLast();
1038 TVector *ptrk_p=(TVector*)trlist->
1040 TVector &ptrk=*ptrk_p;
1041 Int_t last_track=Int_t(ptrk(0));
1042 if (last_track==-1) {
1045 trlist->Add(&trinfo);
1047 // check the track list
1048 Int_t nptracks=trlist->GetEntriesFast();
1050 for (Int_t tr=0;tr<nptracks;tr++) {
1051 TVector *pptrk_p=(TVector*)trlist->At(tr);
1052 TVector &pptrk=*pptrk_p;
1053 trk[tr]=Int_t(pptrk(0));
1054 chtrk[tr]=Int_t(pptrk(1));
1056 } // end if nptracks
1058 } //end loop over clusters
1061 //Int_t nentr2=list->GetEntriesFast();
1062 //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
1063 TTree *fAli=gAlice->TreeK();
1066 if (fAli) file =fAli->GetCurrentFile();
1072 //cout<<"start filling digits \n "<<endl;
1073 // const Float_t zero_supm = 6.;
1074 Int_t nentries=list->GetEntriesFast();
1075 //printf(" \n \n nentries %d \n",nentries);
1076 // start filling the digits
1078 for (Int_t nent=0;nent<nentries;nent++) {
1079 AliMUONlist *address=(AliMUONlist*)list->At(nent);
1080 if (address==0) continue;
1081 Int_t ich=address->fChamber;
1082 Int_t q=address->fSignal;
1083 iChamber=(AliMUONchamber*) (*fChambers)[ich];
1084 AliMUONresponse * response=iChamber->GetResponseModel();
1085 Int_t adcmax= (Int_t) response->MaxAdc();
1086 // add white noise and do zero-suppression and signal truncation
1087 Float_t MeanNoise = gRandom->Gaus(1, 0.2);
1088 Float_t Noise = gRandom->Gaus(0, MeanNoise);
1090 if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise;
1091 if ( q <= zero_supm ) continue;
1092 if ( q > adcmax) q=adcmax;
1093 digits[0]=address->fPadX;
1094 digits[1]=address->fPadY;
1096 digits[3]=address->fPhysics;
1097 digits[4]=address->fHit;
1098 //printf("fSignal, fPhysics fTrack %d %d %d \n",digits[2],digits[3],digits[4]);
1100 TObjArray* trlist=(TObjArray*)address->TrackList();
1101 Int_t nptracks=trlist->GetEntriesFast();
1102 //printf("nptracks, trlist %d %p\n",nptracks,trlist);
1104 // this was changed to accomodate the real number of tracks
1105 if (nptracks > 10) {
1106 cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
1110 printf("Attention - nptracks > 2 %d \n",nptracks);
1111 printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
1113 for (Int_t tr=0;tr<nptracks;tr++) {
1114 TVector *pp_p=(TVector*)trlist->At(tr);
1115 if(!pp_p ) printf("pp_p - %p\n",pp_p);
1117 tracks[tr]=Int_t(pp(0));
1118 charges[tr]=Int_t(pp(1));
1119 //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]);
1120 } //end loop over list of tracks for one pad
1121 // Sort list of tracks according to charge
1123 SortTracks(tracks,charges,nptracks);
1125 if (nptracks < 10 ) {
1126 for (Int_t i=nptracks; i<10; i++) {
1133 MUON->AddDigits(ich,tracks,charges,digits);
1137 //cout<<"I'm out of the loops for digitisation"<<endl;
1138 // gAlice->GetEvent(nev);
1139 gAlice->TreeD()->Fill();
1140 //TTree *TD=gAlice->TreeD();
1142 Stat_t ndig=TD->GetEntries();
1143 cout<<"number of digits "<<ndig<<endl;
1145 for (int i=0;i<10;i++) {
1146 fDch= MUON->DigitsAddress(i);
1147 int ndig=fDch->GetEntriesFast();
1148 printf (" i, ndig %d %d \n",i,ndig);
1151 MUON->ResetDigits();
1154 for(Int_t ii=0;ii<10;++ii) {
1162 } //end loop over cathodes
1165 sprintf(hname,"TreeD%d",nev);
1166 gAlice->TreeD()->Write(hname);
1168 gAlice->TreeD()->Reset();
1170 //Int_t nadr=p_adr->GetEntriesFast();
1171 // printf(" \n \n nadr %d \n",nadr);
1173 // start filling the digits
1175 for (Int_t nent=0;nent<nadr;nent++) {
1176 TVector *pv=(TVector*)p_adr->At(nent);
1182 // gObjectTable->Print();
1186 void AliMUON::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
1189 // Sort the list of tracks contributing to a given digit
1190 // Only the 3 most significant tracks are acctually sorted
1194 // Loop over signals, only 3 times
1199 Int_t idx[3] = {-2,-2,-2};
1200 Int_t jch[3] = {-2,-2,-2};
1201 Int_t jtr[3] = {-2,-2,-2};
1204 if (ntr<3) imax=ntr;
1206 for(i=0;i<imax;i++){
1212 if((i == 1 && j == idx[i-1])
1213 ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
1215 if(charges[j] > qmax) {
1223 jch[i]=charges[jmax];
1224 jtr[i]=tracks[jmax];
1241 void AliMUON::FindClusters(Int_t nev,Int_t last_entry)
1245 // Loop on chambers and on cathode planes
1247 for (Int_t icat=0;icat<2;icat++) {
1248 gAlice->ResetDigits();
1249 gAlice->TreeD()->GetEvent(last_entry+icat); // spurious +1 ...
1250 if (nev < 10) printf("last_entry , icat - %d %d \n",last_entry,icat);
1251 //gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
1253 for (Int_t ich=0;ich<10;ich++) {
1254 AliMUONchamber* iChamber=(AliMUONchamber*) (*fChambers)[ich];
1255 TClonesArray *MUONdigits = this->DigitsAddress(ich);
1256 if (MUONdigits == 0) continue;
1258 // Get ready the current chamber stuff
1260 AliMUONresponse* response = iChamber->GetResponseModel();
1261 AliMUONsegmentation* seg = iChamber->GetSegmentationModel(icat+1);
1262 AliMUONClusterFinder* rec = iChamber->GetReconstructionModel();
1263 // if (icat==1 && (ich==4 || ich==5)) continue;
1264 // printf("icat, ich, seg - %d %d %p\n",icat,ich,seg);
1266 rec->SetSegmentation(seg);
1267 rec->SetResponse(response);
1268 rec->SetDigits(MUONdigits);
1269 rec->SetChamber(ich);
1270 if (nev==0) rec->CalibrateCOG();
1271 // rec->CalibrateCOG();
1272 rec->FindRawClusters();
1274 //printf("Finish FindRawClusters for cathode %d in chamber %d\n",icat,ich);
1277 fRch=RawClustAddress(ich);
1284 //TTree *TR=gAlice->TreeR();
1286 gAlice->TreeR()->Fill();
1288 Stat_t nent=TR->GetEntries();
1289 cout<<"number of entries "<<nent<<endl;
1291 for (int i=0;i<10;i++) {
1292 fRch=RawClustAddress(i);
1293 int nraw=fRch->GetEntriesFast();
1294 printf (" i, nraw %d %d \n",i,nraw);
1302 sprintf(hname,"TreeR%d",nev);
1303 gAlice->TreeR()->Write(hname);
1304 gAlice->TreeR()->Reset();
1306 //gObjectTable->Print();
1310 //______________________________________________________________________________
1311 //_____________________________________________________________________________
1312 void AliMUON::CathodeCorrelation(Int_t nev)
1315 // Correlates the clusters on the two cathode planes and build a list of
1316 // other possible combinations (potential ghosts) - for the moment use the
1317 // criteria of minimum distance between the CoGs of the two correlated
1322 // Loop on chambers and on clusters on the cathode plane with the highest
1323 // number of clusters
1325 static Bool_t first=kTRUE;
1327 AliMUONRawCluster *mRaw1;
1328 AliMUONRawCluster *mRaw2;
1329 AliMUONchamber *iChamber;
1330 AliMUONsegmentation *seg;
1331 TArrayF x1, y1, x2, y2, q1, q2;
1339 // Get pointers to Alice detectors and Digits containers
1340 TTree *TR = gAlice->TreeR();
1341 Int_t nent=(Int_t)TR->GetEntries();
1342 if (nev < 10) printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
1346 Float_t xc2[4],yc2[4];
1347 Float_t xrec2, yrec2;
1348 Float_t xd0, xdif, ydif;
1349 Float_t ysrch,xd,xmax,ymax;
1350 Int_t ilow, iup, iraw1, i;
1353 Float_t xdarray[50];
1358 // Int_t nraw[2], entry,cathode;
1360 for (i=0;i<50;i++) {
1373 // access to the Raw Clusters tree
1374 for (Int_t ich=0;ich<10;ich++) {
1375 iChamber = &(Chamber(ich));
1376 TClonesArray *MUONrawclust = RawClustAddress(ich);
1378 TR->GetEvent(nent-2);
1380 Int_t nrawcl1 = MUONrawclust->GetEntries();
1381 // printf("Found %d raw clusters for cathode 1 in chamber %d \n"
1383 if (!nrawcl1) continue;
1385 seg = iChamber->GetSegmentationModel(1);
1386 // loop over raw clusters of first cathode
1387 for (iraw1=0; iraw1<nrawcl1; iraw1++) {
1388 mRaw1= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw1);
1389 x1[iraw1]=mRaw1->fX;
1390 y1[iraw1]=mRaw1->fY;
1391 q1[iraw1]=(Float_t)mRaw1->fQ; //maybe better fPeakSignal
1392 } // rawclusters cathode 1
1394 // Get information from 2nd cathode
1396 TR->GetEvent(nent-1);
1398 Int_t nrawcl2 = MUONrawclust->GetEntries();
1400 for (iraw1=0; iraw1<nrawcl1; iraw1++) {
1404 //printf("nrawcl2 is zero - idx[0] %d\n",idx[0]);
1406 AddCathCorrel(ich,idx,xc2,yc2);
1412 } // store information from cathode 1 only
1414 // printf("Found %d raw clusters for cathode 2 in chamber %d \n",
1417 for (Int_t iraw2=0; iraw2<nrawcl2; iraw2++) {
1418 mRaw2= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw2);
1419 x2[iraw2]=mRaw2->fX;
1420 y2[iraw2]=mRaw2->fY;
1421 q2[iraw2]=(Float_t)mRaw2->fQ;
1422 } // rawclusters cathode 2
1424 // Initalisation finished
1425 for (iraw1=0; iraw1<nrawcl1; iraw1++) {
1428 seg->GetPadIxy(x1[iraw1],y1[iraw1],ix,iy);
1429 Int_t isec=seg->Sector(ix,iy);
1430 // range to look for ghosts ?!
1432 ymax = seg->Dpy(isec)*7/2;
1433 xmax = seg->Dpx(isec)*7/2;
1435 ymax = seg->Dpy(isec)*13/2;
1436 xmax = seg->Dpx(isec)*3/2;
1438 ysrch=ymax+y1[iraw1];
1440 ilow = AliMUONRawCluster::
1441 BinarySearch(ysrch-2*ymax,y2,0,nrawcl2+1);
1442 iup= AliMUONRawCluster::
1443 BinarySearch(ysrch,y2,ilow,nrawcl2+1);
1444 if (ilow<0 || iup <0 || iup>nrawcl2) continue;
1446 for (Int_t iraw2=ilow; iraw2<=iup; iraw2++) {
1449 xdif=x1[iraw1]-xrec2;
1450 ydif=y1[iraw1]-yrec2;
1451 xd=TMath::Sqrt(xdif*xdif+ydif*ydif);
1455 Sqrt(2*xmax*2*xmax+2*ymax*2*ymax);
1458 Float_t qdif=TMath::Abs(q1[iraw1]-q2[iraw2])/q1[iraw1];
1460 if (x1[iraw1]*xrec2 > 0) {
1462 // printf("q1, q2 qdif % f %f %f \n",q1[iraw1],q2[iraw2],qdif);
1463 // printf("x1, x2 y1 y2 % f %f %f %f \n",x1[iraw1],xrec2,y1[iraw1],yrec2);
1464 //if (qdif <0.3) { //check this number
1467 idx2[counter]=iraw2;
1468 xdarray[counter]=xd;
1469 xarray[counter]=xdif;
1470 yarray[counter]=ydif;
1471 qarray[counter]=qdif;
1476 } // check for same quadrant
1477 } // loop over 2nd cathode range
1482 SortMin(idx2,xdarray,xarray,yarray,qarray,counter);
1483 if (xdarray[0]<seg->Dpx(isec) && xdarray[1]<seg->Dpx(isec)) {
1484 if (qarray[0]>qarray[1]){
1492 if (counter <3) imax=counter;
1495 for (int i=0;i<imax;i++) {
1496 if (idx2[i] >= 0 && idx2[i] < nrawcl2) {
1497 if (xarray[i] > xmax || yarray[i] > 2*ymax)
1504 // add info about the cluster on the 'starting' cathode
1509 //if (idx[0] <0) printf("iraw1 imax idx2[0] idx[0] %d %d %d %d\n",iraw1,imax,idx2[0],idx[0]);
1510 AddCathCorrel(ich,idx,xc2,yc2);
1512 for (Int_t ii=0;ii<counter;ii++) {
1519 for (Int_t iii=0;iii<3;iii++) {
1540 //Int_t nentries=(Int_t)TC->GetEntries();
1541 //cout<<"number entries in tree of correlated clusters "<<nentries<<endl;
1543 static Int_t countev=0;
1546 for (Int_t ii=0;ii<10;ii++) {
1547 fCch= CathCorrelAddress(ii);
1548 Int_t ncor=fCch->GetEntriesFast();
1549 printf (" ii, ncor %d %d \n",ii,ncor);
1550 if (ncor>=2) countch++;
1555 sprintf(hname,"TreeC%d",nev);
1561 if (countch==10) countev++;
1562 printf("countev - %d\n",countev);
1564 // gObjectTable->Print();
1570 //_____________________________________________________________________________
1572 void AliMUON::MakeTreeC(Option_t *option)
1574 char *C = strstr(option,"C");
1575 if (C && !fTreeC) fTreeC = new TTree("TC","CathodeCorrelation");
1577 // Create a branch for correlation
1579 const Int_t buffersize = 4000;
1580 char branchname[30];
1582 // one branch for correlation per chamber
1583 for (int i=0; i<10 ;i++) {
1584 sprintf(branchname,"%sCorrelation%d",GetName(),i+1);
1586 if (fCathCorrel && fTreeC) {
1587 TreeC()->Branch(branchname,&((*fCathCorrel)[i]), buffersize);
1588 printf("Making Branch %s for correlation in chamber %d\n",branchname,i+1);
1593 //_____________________________________________________________________________
1594 void AliMUON::GetTreeC(Int_t event)
1597 // set the branch address
1599 char branchname[30];
1606 sprintf(treeName,"TreeC%d",event);
1607 fTreeC = (TTree*)gDirectory->Get(treeName);
1612 for (int i=0; i<10; i++) {
1613 sprintf(branchname,"%sCorrelation%d",GetName(),i+1);
1615 branch = fTreeC->GetBranch(branchname);
1616 if (branch) branch->SetAddress(&((*fCathCorrel)[i]));
1620 printf("ERROR: cannot find CathodeCorrelation Tree for event:%d\n",event);
1623 // gObjectTable->Print();
1628 void AliMUON::Streamer(TBuffer &R__b)
1630 // Stream an object of class AliMUON.
1631 AliMUONchamber *iChamber;
1632 AliMUONsegmentation *segmentation;
1633 AliMUONresponse *response;
1634 TClonesArray *digitsaddress;
1635 TClonesArray *rawcladdress;
1636 TClonesArray *corcladdress;
1637 // TObjArray *clustaddress;
1639 if (R__b.IsReading()) {
1640 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1641 AliDetector::Streamer(R__b);
1643 R__b >> fClusters; // diff
1645 R__b >> fRawClusters;
1646 R__b >> fCathCorrel;
1647 R__b.ReadArray(fNdch);
1648 R__b.ReadArray(fNrawch);
1649 R__b.ReadArray(fNcorch);
1650 // R__b >> fRecClusters;
1663 // Stream chamber related information
1664 for (Int_t i =0; i<10; i++) {
1665 iChamber=(AliMUONchamber*) (*fChambers)[i];
1666 iChamber->Streamer(R__b);
1667 if (iChamber->Nsec()==1) {
1668 segmentation=iChamber->GetSegmentationModel(1);
1669 segmentation->Streamer(R__b);
1671 segmentation=iChamber->GetSegmentationModel(1);
1672 segmentation->Streamer(R__b);
1673 segmentation=iChamber->GetSegmentationModel(2);
1674 segmentation->Streamer(R__b);
1676 response=iChamber->GetResponseModel();
1677 response->Streamer(R__b);
1678 digitsaddress=(TClonesArray*) (*fDchambers)[i];
1679 digitsaddress->Streamer(R__b);
1680 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1681 rawcladdress->Streamer(R__b);
1682 corcladdress=(TClonesArray*) (*fCathCorrel)[i];
1683 corcladdress->Streamer(R__b);
1687 R__b.WriteVersion(AliMUON::IsA());
1688 AliDetector::Streamer(R__b);
1690 R__b << fClusters; // diff
1692 R__b << fRawClusters;
1693 R__b << fCathCorrel;
1694 R__b.WriteArray(fNdch, 10);
1695 R__b.WriteArray(fNrawch, 10);
1696 R__b.WriteArray(fNcorch, 10);
1697 // R__b << fRecClusters;
1710 // Stream chamber related information
1712 for (Int_t i =0; i<20; i++) {
1713 clustaddress=(TObjArray*) (*fRecClusters)[i];
1714 clustaddress->Streamer(R__b);
1717 for (Int_t i =0; i<10; i++) {
1718 iChamber=(AliMUONchamber*) (*fChambers)[i];
1719 iChamber->Streamer(R__b);
1720 if (iChamber->Nsec()==1) {
1721 segmentation=iChamber->GetSegmentationModel(1);
1722 segmentation->Streamer(R__b);
1724 segmentation=iChamber->GetSegmentationModel(1);
1725 segmentation->Streamer(R__b);
1726 segmentation=iChamber->GetSegmentationModel(2);
1727 segmentation->Streamer(R__b);
1729 response=iChamber->GetResponseModel();
1730 response->Streamer(R__b);
1731 digitsaddress=(TClonesArray*) (*fDchambers)[i];
1732 digitsaddress->Streamer(R__b);
1733 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1734 rawcladdress->Streamer(R__b);
1735 corcladdress=(TClonesArray*) (*fCathCorrel)[i];
1736 corcladdress->Streamer(R__b);
1740 AliMUONcluster* AliMUON::FirstPad(AliMUONhit* hit, TClonesArray *clusters)
1743 // Initialise the pad iterator
1744 // Return the address of the first padhit for hit
1745 TClonesArray *theClusters = clusters;
1746 Int_t nclust = theClusters->GetEntriesFast();
1747 if (nclust && hit->fPHlast > 0) {
1748 sMaxIterPad=hit->fPHlast;
1749 sCurIterPad=hit->fPHfirst;
1750 return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1);
1756 AliMUONcluster* AliMUON::NextPad(TClonesArray *clusters)
1759 if (sCurIterPad <= sMaxIterPad) {
1760 return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1);
1766 //////////////////////////// modifs perso ///////////////
1768 static TTree *ntuple_global;
1769 static TFile *hfile_global;
1771 // variables of the tracking ntuple
1773 Int_t ievr; // number of event
1774 Int_t ntrackr; // number of tracks per event
1775 Int_t istatr[500]; // 1 = good muon, 2 = ghost, 0 = something else
1776 Int_t isignr[500]; // sign of the track
1777 Float_t pxr[500]; // x momentum of the reconstructed track
1778 Float_t pyr[500]; // y momentum of the reconstructed track
1779 Float_t pzr[500]; // z momentum of the reconstructed track
1780 Float_t zvr[500]; // z vertex
1781 Float_t chi2r[500]; // chi2 of the fit of the track with the field map
1782 Float_t pxv[500]; // x momentum at vertex
1783 Float_t pyv[500]; // y momentum at vertex
1784 Float_t pzv[500]; // z momentum at vertex
1787 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
1789 TClonesArray *MUONrawclust = RawClustAddress(ichamber);
1791 TTree *TR = gAlice->TreeR();
1792 Int_t nent=(Int_t)TR->GetEntries();
1793 TR->GetEvent(nent-2+icathod-1);
1794 //TR->GetEvent(icathod);
1795 //Int_t nrawcl = (Int_t)MUONrawclust->GetEntriesFast();
1797 AliMUONRawCluster * mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(icluster);
1798 //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
1803 void 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)
1806 // open kine and hits tree of background file for reconstruction of geant hits
1807 // call tracking fortran program
1808 static Bool_t first=kTRUE;
1810 char *Add = strstr(option,"Add");
1812 if (Add ) { // only in case of background with geant hits
1815 cout<<"filename "<<fFileName<<endl;
1816 File=new TFile(fFileName);
1817 cout<<"I have opened "<<fFileName<<" file "<<endl;
1818 fHits2 = new TClonesArray("AliMUONhit",1000 );
1819 fParticles2 = new TClonesArray("GParticle",1000);
1823 if(fHits2) fHits2->Clear();
1824 if(fParticles2) fParticles2->Clear();
1825 if(TrH1) delete TrH1;
1829 // Get Hits Tree header from file
1831 sprintf(treeName,"TreeH%d",bgd_ev);
1832 TrH1 = (TTree*)gDirectory->Get(treeName);
1834 printf("ERROR: cannot find Hits Tree for event:%d\n",bgd_ev);
1836 // set branch addresses
1838 char branchname[30];
1839 sprintf(branchname,"%s",GetName());
1840 if (TrH1 && fHits2) {
1841 branch = TrH1->GetBranch(branchname);
1842 if (branch) branch->SetAddress(&fHits2);
1845 // get the Kine tree
1846 sprintf(treeName,"TreeK%d",bgd_ev);
1847 TK1 = (TTree*)gDirectory->Get(treeName);
1849 printf("ERROR: cannot find Kine Tree for event:%d\n",bgd_ev);
1851 // set branch addresses
1853 TK1->SetBranchAddress("Particles", &fParticles2);
1856 // get back to the first file
1857 TTree *TK = gAlice->TreeK();
1859 if (TK) file1 = TK->GetCurrentFile();
1864 // call tracking fortran program
1865 reconstmuon(ifit,idebug,nev,idres,ireadgeant);
1869 void AliMUON::Init(Double_t &seff, Double_t &sb0, Double_t &sbl3)
1872 // introduce in fortran program somme parameters and cuts for tracking
1873 // create output file "reconst.root" (histos + ntuple)
1874 cutpxz(fSPxzCut); // Pxz cut (GeV/c) to begin the track finding
1875 sigmacut(fSSigmaCut); // Number of sigmas delimiting the searching areas
1876 xpreci(fSXPrec); // Chamber precision in X (cm)
1877 ypreci(fSYPrec); // Chamber precision in Y (cm)
1878 reco_init(seff,sb0,sbl3);
1881 void AliMUON::FinishEvent()
1883 TTree *TK = gAlice->TreeK();
1885 if (TK) file1 = TK->GetCurrentFile();
1889 void AliMUON::Close()
1892 // write histos and ntuple to "reconst.root" file
1896 void chfill(Int_t &id, Float_t &x, Float_t &y, Float_t &w)
1899 // fill histo like hfill in fortran
1901 sprintf(name,"h%d",id);
1902 TH1F *h1 = (TH1F*) gDirectory->Get(name);
1906 void chfill2(Int_t &id, Float_t &x, Float_t &y, Float_t &w)
1909 // fill histo like hfill2 in fortran
1911 sprintf(name,"h%d",id);
1912 TH2F *h2 = (TH2F*) gDirectory->Get(name);
1916 void chf1(Int_t &id, Float_t &x, Float_t &w)
1919 // fill histo like hf1 in fortran
1921 sprintf(name,"h%d",id);
1922 TH1F *h1 = (TH1F*) gDirectory->Get(name);
1929 // Create an output file ("reconst.root")
1930 // Create some histograms and an ntuple
1932 hfile_global = new TFile("reconst.root","RECREATE","Ntuple - reconstruction");
1934 ntuple_global = new TTree("ntuple","Reconst ntuple");
1935 ntuple_global->Branch("ievr",&ntuple_st.ievr,"ievr/I");
1936 ntuple_global->Branch("ntrackr",&ntuple_st.ntrackr,"ntrackr/I");
1937 ntuple_global->Branch("istatr",&ntuple_st.istatr[0],"istatr[500]/I");
1938 ntuple_global->Branch("isignr",&ntuple_st.isignr[0],"isignr[500]/I");
1939 ntuple_global->Branch("pxr",&ntuple_st.pxr[0],"pxr[500]/F");
1940 ntuple_global->Branch("pyr",&ntuple_st.pyr[0],"pyr[500]/F");
1941 ntuple_global->Branch("pzr",&ntuple_st.pzr[0],"pzr[500]/F");
1942 ntuple_global->Branch("zvr",&ntuple_st.zvr[0],"zvr[500]/F");
1943 ntuple_global->Branch("chi2r",&ntuple_st.chi2r[0],"chi2r[500]/F");
1944 ntuple_global->Branch("pxv",&ntuple_st.pxv[0],"pxv[500]/F");
1945 ntuple_global->Branch("pyv",&ntuple_st.pyv[0],"pyv[500]/F");
1946 ntuple_global->Branch("pzv",&ntuple_st.pzv[0],"pzv[500]/F");
1950 new TH1F("h100","particule id du hit geant",20,0.,20.);
1951 new TH1F("h101","position en x du hit geant",100,-200.,200.);
1952 new TH1F("h102","position en y du hit geant",100,-200.,200.);
1953 new TH1F("h103","chambre de tracking concernee",15,0.,14.);
1954 new TH1F("h104","moment ptot du hit geant",50,0.,100.);
1955 new TH1F("h105","px au vertex",50,0.,20.);
1956 new TH1F("h106","py au vertex",50,0.,20.);
1957 new TH1F("h107","pz au vertex",50,0.,20.);
1958 new TH1F("h108","position zv",50,-15.,15.);
1959 new TH1F("h109","position en x du hit reconstruit",100,-300.,300.);
1960 new TH1F("h110","position en y du hit reconstruit",100,-300.,300.);
1961 new TH1F("h111","delta x ",100,-0.4,0.4);
1962 new TH1F("h112","delta y ",100,-0.4,0.4);
1966 for (int i=0;i<10;i++) {
1967 sprintf(hname,"deltax%d",i);
1968 sprintf(hname1,"h12%d",i);
1969 new TH1F(hname1,hname ,100,-0.4,0.4);
1970 sprintf(hname,"deltay%d",i);
1971 sprintf(hname1,"h13%d",i);
1972 new TH1F(hname1,hname ,100,-0.4,0.4);
1974 new TH2F("h2000","VAR X st. 5",30,3.0,183.0,100,0.,25.);
1975 new TH2F("h2001","VAR Y st. 5",30,3.0,183.0,100,0.,25.);
1977 new TH2F("h2500","P vs X HHIT",30,3.0,183.0,200,0.,200.);
1978 new TH2F("h2501","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.);
1979 new TH2F("h2502","P vs X EPH2 st. 5",30,3.0,183.0,100,0.,0.000005);
1980 new TH2F("h2503","P vs X EAL2 st. 5",30,3.0,183.0,100,0.,0.01);
1981 //new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,1.5);
1982 new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,0.1);
1983 new TH2F("h2505","P vs X EYM2 st. 5",30,3.0,183.0,100,0.,30.);
1985 new TH2F("h2507","P vs X EPH st. 5",30,3.0,183.0,100,0.,0.003);
1986 new TH2F("h2508","P vs X EAL st. 5",30,3.0,183.0,100,0.,0.3);
1987 //new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,1.5);
1988 new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,0.4);
1989 new TH2F("h2510","P vs X EYM st. 5",30,3.0,183.0,100,0.,30.);
1991 new TH2F("h2511","P vs X EPH cut st. 5",30,3.0,183.0,100,0.,0.01);
1992 new TH2F("h2512","P vs X EAL cut st. 5",30,3.0,183.0,100,0.,0.3);
1993 //new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,1.5);
1994 new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,0.4);
1995 new TH2F("h2514","P vs X EYM cut st. 5",30,3.0,183.0,100,0.,30.);
1997 new TH2F("h2400","P vs X HHIT",30,3.0,183.0,200,0.,200.);
1998 new TH2F("h2401","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.);
1999 new TH2F("h2402","P vs X EPH2 st. 4",30,3.0,183.0,100,0.,0.000005);
2000 new TH2F("h2403","P vs X EAL2 st. 4",30,3.0,183.0,100,0.,0.05);
2001 //new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,1.5);
2002 new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,0.1);
2003 new TH2F("h2405","P vs X EYM2 st. 4",30,3.0,183.0,100,0.,30.);
2005 new TH2F("h2407","P vs X EPH st. 4",30,3.0,183.0,100,0.,0.003);
2006 new TH2F("h2408","P vs X EAL st. 4",30,3.0,183.0,100,0.,0.3);
2007 //new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,1.5);
2008 new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,0.1);
2009 new TH2F("h2410","P vs X EYM st. 4",30,3.0,183.0,100,0.,30.);
2011 new TH2F("h2411","P vs X EPH cut st. 4",30,3.0,183.0,100,0.,0.01);
2012 new TH2F("h2412","P vs X EAL cut st. 4",30,3.0,183.0,100,0.,0.3);
2013 //new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,1.5);
2014 new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,0.1);
2015 new TH2F("h2414","P vs X EYM cut st. 4",30,3.0,183.0,100,0.,30.);
2017 new TH1F("h2301","P2",30,3.0,183.0);
2018 new TH2F("h2302","P2 vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006);
2019 new TH2F("h2303","P2 vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.0005);
2020 //new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5);
2021 new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.);
2022 new TH2F("h2305","P2 vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.);
2024 new TH2F("h2307","P vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006);
2025 new TH2F("h2308","P vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.005);
2026 //new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5);
2027 new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.);
2028 new TH2F("h2310","P vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.);
2030 new TH2F("h2311","P vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06);
2031 new TH2F("h2312","P vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05);
2032 //new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5);
2033 new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.);
2034 new TH2F("h2314","P vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.);
2036 new TH2F("h2315","P2 vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06);
2037 new TH2F("h2316","P2 vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05);
2038 //new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5);
2039 new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.);
2040 new TH2F("h2318","P2 vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.);
2043 new TH1F("h2201","P2",30,3.0,183.0);
2044 new TH2F("h2202","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2045 new TH2F("h2203","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2046 //new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2047 new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2048 new TH2F("h2205","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.);
2050 new TH2F("h2207","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2051 new TH2F("h2208","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2052 //new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2053 new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2054 new TH2F("h2210","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.);
2056 new TH2F("h2211","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05);
2057 new TH2F("h2212","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2058 //new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2059 new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2060 new TH2F("h2214","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.);
2062 new TH2F("h2215","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05);
2063 new TH2F("h2216","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2064 //new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2065 new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2066 new TH2F("h2218","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.);
2069 new TH2F("h2102","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2070 new TH2F("h2103","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2071 //new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2072 new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2073 new TH2F("h2105","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.);
2075 new TH2F("h2107","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2076 new TH2F("h2108","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2077 //new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2078 new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2079 new TH2F("h2110","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.);
2081 new TH2F("h2111","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1);
2082 new TH2F("h2112","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2083 //new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2084 new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2085 new TH2F("h2114","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.);
2087 new TH2F("h2115","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1);
2088 new TH2F("h2116","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2089 //new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2090 new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2091 new TH2F("h2118","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.);
2094 new TH1F("h2701","P2 fit 2",30,3.0,183.0);
2095 new TH2F("h2702","P2 vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006);
2096 new TH2F("h2703","P2 vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005);
2097 // new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2098 new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.);
2099 new TH2F("h2705","P2 vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.);
2101 new TH2F("h2707","P vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006);
2102 new TH2F("h2708","P vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005);
2103 //new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2104 new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.);
2105 new TH2F("h2710","P vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.);
2107 new TH2F("h2711","P vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07);
2108 new TH2F("h2712","P vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2);
2109 //new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2110 new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.);
2111 new TH2F("h2714","P vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.);
2113 new TH2F("h2715","P2 vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07);
2114 new TH2F("h2716","P2 vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2);
2115 //new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2116 new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.);
2117 new TH2F("h2718","P2 vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.);
2120 new TH1F("h2801","P2 fit 1",30,3.0,183.0);
2121 new TH2F("h2802","P2 vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006);
2122 new TH2F("h2803","P2 vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005);
2123 //new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2124 new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.);
2125 new TH2F("h2805","P2 vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.);
2127 new TH2F("h2807","P vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006);
2128 new TH2F("h2808","P vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005);
2129 //new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2130 new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.);
2131 new TH2F("h2810","P vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.);
2133 new TH2F("h2811","P vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05);
2134 new TH2F("h2812","P vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2);
2135 //new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2136 new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.);
2137 new TH2F("h2814","P vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.);
2139 new TH2F("h2815","P2 vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05);
2140 new TH2F("h2816","P2 vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2);
2141 //new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2142 new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.);
2143 new TH2F("h2818","P2 vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.);
2146 new TH1F("h500","Acceptance en H st. 4",500,0.,500.);
2147 new TH1F("h600","Acceptance en H st. 5",500,0.,500.);
2148 new TH1F("h700","X vertex track found",200,-10.,10.);
2149 new TH1F("h701","Y vertex track found",200,-10.,10.);
2150 new TH1F("h800","Rap. muon gen.",100,0.,5.);
2151 new TH1F("h801","Rap. muon gen. recons.",100,0.,5.);
2152 new TH1F("h802","Rap. muon gen. ghost ",100,0.,5.);
2153 new TH1F("h900","Pt muon gen.",100,0.,20.);
2154 new TH1F("h901","Pt muon gen. recons.",100,0.,20.);
2155 new TH1F("h902","Pt muon gen. ghost",100,0.,20.);
2156 new TH1F("h910","phi muon gen.",100,-10.,10.);
2157 new TH1F("h911","phi muon gen. recons.",100,-10.,10.);
2158 new TH1F("h912","phi muon gen. ghost",100,-10.,10.);
2159 new TH2F("h1001","Y VS X hit st. 1",300,-300.,300.,300,-300.,300.);
2160 new TH2F("h1002","Y VS X hit st. 2",300,-300.,300.,300,-300.,300.);
2161 new TH2F("h1003","Y VS X hit st. 3",300,-300.,300.,300,-300.,300.);
2162 new TH2F("h1004","Y VS X hit st. 4",300,-300.,300.,300,-300.,300.);
2163 new TH2F("h1005","Y VS X hit st. 5",300,-300.,300.,300,-300.,300.);
2164 // Histos variance dans 4
2165 new TH2F("h11","VAR X st. 4",30,3.0,183.0,100,0.,2.);
2166 new TH2F("h12","VAR Y st. 4",30,3.0,183.0,100,0.,600.);
2167 new TH2F("h13","VAR PHI st. 4",30,3.0,183.0,100,0.,0.0001);
2168 new TH2F("h14","VAR ALM st. 4",30,3.0,183.0,100,0.,0.05);
2169 new TH1F("h15","P",30,3.0,183.0);
2170 new TH1F("h411","VAR X st. 4",100,-1.42,1.42);
2171 new TH1F("h412","VAR Y st. 4",100,-25.,25.);
2172 new TH1F("h413","VAR PHI st. 4",100,-0.01,0.01);
2173 new TH1F("h414","VAR ALM st. 4",100,-0.23,0.23);
2175 new TH2F("h211","histo2-VAR X st. 4",30,3.0,183.0,100,0.,2.);
2176 new TH2F("h212","histo2-VAR Y st. 4",30,3.0,183.0,100,0.,600.);
2177 new TH1F("h213","histo2-VAR X st. 4",100,-1.42,1.42);
2178 new TH1F("h214","histo2-VAR Y st. 4",100,-25.,25.);
2179 new TH1F("h215","histo2-P",30,3.0,183.0);
2181 // Histos variance dans 2
2182 new TH2F("h21","VAR X st. 2",30,3.0,183.0,100,0.,3.);
2183 new TH2F("h22","VAR Y st. 2",30,3.0,183.0,100,0.,7.);
2184 new TH2F("h23","VAR PHI st. 2",30,3.0,183.0,100,0.,0.006);
2185 new TH2F("h24","VAR ALM st. 2",30,3.0,183.0,100,0.,0.005);
2186 new TH1F("h25","P",30,3.0,183.0);
2187 new TH1F("h421","VAR X st. 2",100,-1.72,1.72);
2188 new TH1F("h422","VAR Y st. 2",100,-2.7,2.7);
2189 new TH1F("h423","VAR PHI st. 2",100,-0.08,0.08);
2190 new TH1F("h424","VAR ALM st. 2",100,-0.072,0.072);
2192 new TH2F("h221","histo2-VAR X st. 2",30,3.0,183.0,100,0.,3.);
2193 new TH2F("h222","histo2-VAR Y st. 2",30,3.0,183.0,100,0.,7.);
2194 new TH1F("h223","histo2-VAR X st. 2",100,-1.72,1.72);
2195 new TH1F("h224","histo2-VAR Y st. 2",100,-2.7,2.7);
2196 new TH1F("h225","histo2-P",30,3.0,183.0);
2198 // Histos variance dans 1
2199 new TH2F("h31","VAR X st. 1",30,3.0,183.0,100,0.,2.);
2200 new TH2F("h32","VAR Y st. 1",30,3.0,183.0,100,0.,0.5);
2201 new TH2F("h33","VAR PHI st. 1",30,3.0,183.0,100,0.,0.006);
2202 new TH2F("h34","VAR ALM st. 1",30,3.0,183.0,100,0.,0.005);
2203 new TH1F("h35","P",30,3.0,183.0);
2204 new TH1F("h431","VAR X st. 1",100,-1.42,1.42);
2205 new TH1F("h432","VAR Y st. 1",100,-0.72,0.72);
2206 new TH1F("h433","VAR PHI st. 1",100,-0.08,0.08);
2207 new TH1F("h434","VAR ALM st. 1",100,-0.072,0.072);
2208 // Histos variance dans 1
2209 new TH2F("h41","VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.);
2210 new TH2F("h42","VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.);
2211 new TH2F("h43","VAR PHI st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005);
2212 new TH2F("h44","VAR ALM st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005);
2213 new TH1F("h45","P",30,3.0,183.0);
2214 new TH1F("h441","VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.);
2215 new TH1F("h442","VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5);
2216 new TH1F("h443","VAR PHI st. 1 fit 5,4,3,2,V",100,-0.072,0.072);
2217 new TH1F("h444","VAR ALM st. 1 fit 5,4,3,2,V",100,-0.072,0.072);
2219 new TH2F("h241","histo2-VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.);
2220 new TH2F("h242","histo2-VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.);
2221 new TH1F("h243","histo2-VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.);
2222 new TH1F("h244","histo2-VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5);
2223 new TH1F("h245","histo2-P",30,3.0,183.0);
2225 // Histos variance dans 2
2226 new TH2F("h51","VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5);
2227 new TH2F("h52","VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.);
2228 new TH2F("h53","VAR PHI st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.005);
2229 new TH2F("h54","VAR ALM st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.01);
2230 new TH1F("h55","P",30,3.0,183.0);
2231 new TH1F("h451","VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72);
2232 new TH1F("h452","VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42);
2233 new TH1F("h453","VAR PHI st. 2 fit 5,4,3,1,V",100,-0.072,0.072);
2234 new TH1F("h454","VAR ALM st. 2 fit 5,4,3,1,V",100,-0.1,0.1);
2235 new TH1F("h999","PTOT",30,3.0,183.0);
2237 new TH2F("h251","histo2-VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5);
2238 new TH2F("h252","histo2-VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.);
2239 new TH1F("h253","histo2-VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72);
2240 new TH1F("h254","histo2-VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42);
2241 new TH1F("h255","histo2-P",30,3.0,183.0);
2242 // Histos variance dans 3
2243 new TH2F("h61","VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.);
2244 new TH2F("h62","VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.);
2245 new TH2F("h63","VAR PHI st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006);
2246 new TH2F("h64","VAR ALM st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006);
2247 new TH1F("h65","P",30,3.0,183.0);
2248 new TH1F("h461","VAR X st. 3 fit 4,5,V",100,-2.25,2.25);
2249 new TH1F("h462","VAR Y st. 3 fit 4,5,V",100,-1.42,1.42);
2250 new TH1F("h463","VAR PHI st. 3 fit 4,5,V",100,-0.024,0.024);
2251 new TH1F("h464","VAR ALM st. 3 fit 4,5,V",100,-0.024,0.024);
2253 new TH2F("h261","histo2-VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.);
2254 new TH2F("h262","histo2-VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.);
2255 new TH1F("h263","histo2-VAR X st. 3 fit 4,5,V",100,-2.25,2.25);
2256 new TH1F("h264","histo2-VAR Y st. 3 fit 4,5,V",100,-1.42,1.42);
2257 new TH1F("h265","Phisto2-",30,3.0,183.0);
2258 // Histos dx,dy distribution between chambers inside stations
2259 new TH1F("h71","DX in st. ID-70",100,-5.,5.);
2260 new TH1F("h81","DY in st. ID-80",100,-5.,5.);
2261 new TH1F("h72","DX in st. ID-70",100,-5.,5.);
2262 new TH1F("h82","DY in st. ID-80",100,-5.,5.);
2263 new TH1F("h73","DX in st. ID-70",100,-5.,5.);
2264 new TH1F("h83","DY in st. ID-80",100,-5.,5.);
2265 new TH1F("h74","DX in st. ID-70",100,-5.,5.);
2266 new TH1F("h84","DY in st. ID-80",100,-5.,5.);
2267 new TH1F("h75","DX in st. ID-70",100,-5.,5.);
2268 new TH1F("h85","DY in st. ID-80",100,-5.,5.);
2271 void 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)
2275 ntuple_st.ievr = ievr;
2276 ntuple_st.ntrackr = ntrackr;
2277 for (Int_t i=0; i<500; i++) {
2278 ntuple_st.istatr[i] = istatr[i];
2279 ntuple_st.isignr[i] = isignr[i];
2280 ntuple_st.pxr[i] = pxr[i];
2281 ntuple_st.pyr[i] = pyr[i];
2282 ntuple_st.pzr[i] = pzr[i];
2283 ntuple_st.zvr[i] = zvr[i];
2284 ntuple_st.chi2r[i] = chi2r[i];
2285 ntuple_st.pxv[i] = pxv[i];
2286 ntuple_st.pyv[i] = pyv[i];
2287 ntuple_st.pzv[i] = pzv[i];
2289 ntuple_global->Fill();
2295 // write histos and ntuple to "reconst.root" file
2296 hfile_global->Write();
2299 void 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)
2302 // introduce aliroot variables in fortran common
2303 // tracking study from geant hits
2306 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
2308 // TTree *TK = gAlice->TreeK();
2309 TTree *TH = gAlice->TreeH();
2310 Int_t ntracks = (Int_t)TH->GetEntries();
2311 cout<<"ntrack="<<ntracks<<endl;
2320 for (Int_t track=0; track<ntracks;track++) {
2321 gAlice->ResetHits();
2322 TH->GetEvent(track);
2328 for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
2330 mHit=(AliMUONhit*)MUON->NextHit())
2332 if (maxidg<=20000) {
2334 if (mHit->fChamber > 10) continue;
2335 TClonesArray *fPartArray = gAlice->Particles();
2337 Int_t ftrack = mHit->fTrack;
2338 Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
2340 if (id==kMuonPlus||id==kMuonMinus) {
2342 // inversion de x et y car le champ est inverse dans le programme tracking
2345 xgeant[maxidg] = mHit->fY; // x-pos of hit
2346 ygeant[maxidg] = mHit->fX; // y-pos of hit
2347 clsize1[maxidg] = 0; // cluster size on 1-st cathode
2348 clsize2[maxidg] = 0; // cluster size on 2-nd cathode
2349 cx[maxidg] = mHit->fCyHit; // Px/P of hit
2350 cy[maxidg] = mHit->fCxHit; // Py/P of hit
2351 cz[maxidg] = mHit->fCzHit; // Pz/P of hit
2352 izch[maxidg] = mHit->fChamber;
2354 Int_t pdgtype = Int_t(mHit->fParticle); // particle number
2355 itypg[maxidg] = gMC->IdFromPDG(pdgtype);
2358 if (id==kMuonPlus) itypg[maxidg] = 5;
2359 else itypg[maxidg] = 6;
2363 //printf("ich, itypg[maxidg] %d %d\n",izch[maxidg],itypg[maxidg]);
2365 ptotg[maxidg] = mHit->fPTot; // P of hit
2367 Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
2368 Float_t thet = Part->Theta();
2369 thet = thet*180./3.1416;
2371 //cout<<"chambre "<<izch[maxidg]<<" ptot="<<ptotg[maxidg]<<" theta="<<thet<<" phi="<<mHit->fPhi<<" z="<<zz<<endl;
2373 Int_t iparent = Part->GetFirstMother();
2377 ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother();
2385 //printf("iparent - %d\n",iparent);
2386 Int_t id1 = ftrack; // numero de la particule generee au vertex
2387 Int_t idum = track+1;
2388 Int_t id2 = ((TParticle*) fPartArray->UncheckedAt(iparent))->GetPdgCode();
2390 if (id2==443) id2=114;
2396 //printf("id2 %d\n",id2);
2397 idg[maxidg] = 30000*id1+10000*idum+id2;
2399 pvert1g[maxidg] = Part->Py(); // Px vertex
2400 pvert2g[maxidg] = Part->Px(); // Py vertex
2401 pvert3g[maxidg] = Part->Pz(); // Pz vertex
2402 zvertg[maxidg] = Part->Vz(); // z vertex
2404 // cout<<"x="<<xgeant[maxidg]<<endl;
2405 //cout<<"y="<<ygeant[maxidg]<<endl;
2406 //cout<<"typ="<<itypg[maxidg]<<endl;
2414 } // track loop first file
2416 if (TrH1 && fHits2 ) { // if background file
2417 ntracks =(Int_t)TrH1->GetEntries();
2418 printf("Trackf_read - 2-nd file - ntracks %d\n",ntracks);
2421 for (Int_t track=0; track<ntracks; track++) {
2423 if (fHits2) fHits2->Clear();
2424 TrH1->GetEvent(track);
2427 for (int i=0;i<fHits2->GetEntriesFast();i++)
2429 AliMUONhit *mHit=(AliMUONhit*) (*fHits2)[i];
2431 if (mHit->fChamber > 10) continue;
2433 if (maxidg<=20000) {
2435 // inversion de x et y car le champ est inverse dans le programme tracking !!!!
2436 xtrg[maxidg] = 0; // only for reconstructed point
2437 ytrg[maxidg] = 0; // only for reconstructed point
2438 xgeant[maxidg] = mHit->fY; // x-pos of hit
2439 ygeant[maxidg] = mHit->fX; // y-pos of hit
2440 clsize1[maxidg] = 0; // cluster size on 1-st cathode
2441 clsize2[maxidg] = 0; // cluster size on 2-nd cathode
2442 cx[maxidg] = mHit->fCyHit; // Px/P of hit
2443 cy[maxidg] = mHit->fCxHit; // Py/P of hit
2444 cz[maxidg] = mHit->fCzHit; // Pz/P of hit
2445 izch[maxidg] = mHit->fChamber; // chamber number
2446 ptotg[maxidg] = mHit->fPTot; // P of hit
2448 Int_t ftrack = mHit->fTrack;
2449 Int_t id1 = ftrack; // track number
2450 Int_t idum = track+1;
2452 TClonesArray *fPartArray = fParticles2;
2454 Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
2455 if (id==kMuonPlus||id==kMuonMinus) {
2456 if (id==kMuonPlus) itypg[maxidg] = 5;
2457 else itypg[maxidg] = 6;
2458 } else itypg[maxidg]=0;
2460 Int_t id2=0; // set parent to 0 for background !!
2461 idg[maxidg] = 30000*id1+10000*idum+id2;
2463 pvert1g[maxidg] = Part->Py(); // Px vertex
2464 pvert2g[maxidg] = Part->Px(); // Py vertex
2465 pvert3g[maxidg] = Part->Pz(); // Pz vertex
2466 zvertg[maxidg] = Part->Vz(); // z vertex
2470 } // check limits (maxidg)
2477 cout<<"nhittot1="<<nhittot1<<endl;
2479 static Int_t nbres=0;
2480 if (nres>=19) nbres++;
2481 printf("nres, nbres %d %d \n",nres,nbres);
2489 void 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)
2493 // introduce aliroot variables in fortran common
2494 // tracking study from reconstructed points
2496 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
2498 cout<<"numero de l'evenement "<<nev<<endl;
2500 MUON->GetTreeC(nev);
2501 TTree *TC=MUON->TreeC();
2507 static Int_t nuncor=0;
2508 static Int_t nbadcor=0;
2509 AliMUONRawCluster * mRaw;
2510 AliMUONRawCluster * mRaw1;
2511 TTree *TH = gAlice->TreeH();
2516 for (Int_t ich=0;ich<10;ich++) {
2517 TClonesArray *MUONcorrel = MUON->CathCorrelAddress(ich);
2518 MUON->ResetCorrelation();
2520 Int_t ncor = (Int_t)MUONcorrel->GetEntries();
2521 if (ncor>=2) nncor++;
2522 if (!ncor) continue;
2524 // Loop over correlated clusters
2525 for (Int_t icor=0;icor<ncor;icor++) {
2526 AliMUONcorrelation * mCor = (AliMUONcorrelation*)MUONcorrel->UncheckedAt(icor);
2528 Int_t flag=0; // = 1 if no information in the second cathode
2529 Int_t index = mCor->fCorrelIndex[0]; // for the second cathode
2531 Int_t index1 = mCor->fCorrelIndex[3]; // for the 1-st cathode
2532 mRaw1 = MUON->RawCluster(ich,1,index1);
2533 mult1=mRaw1->fMultiplicity;
2534 mRaw = MUON->RawCluster(ich,2,index);
2535 mult2=mRaw->fMultiplicity;
2537 index = mCor->fCorrelIndex[3];
2538 mRaw = MUON->RawCluster(ich,1,index);
2539 mult1=mRaw->fMultiplicity;
2544 if (!mRaw) continue;
2546 Int_t ftrack1 = mRaw->fTracks[1]; // qui doit etre le meme pour
2547 // la cathode 1 et 2
2548 ihit= mRaw->fTracks[0];
2549 //printf("icor, ftrack1 ihit %d %d %d\n",icor,ftrack1,ihit);
2551 if (mRaw->fClusterType == 0 ) {
2553 if (maxidg<=20000) {
2555 xtrg[maxidg] = (Double_t) mCor->fY[3];
2556 ytrg[maxidg] = (Double_t) mCor->fX[0];
2557 Int_t index1 = mCor->fCorrelIndex[3];
2558 mRaw1 = MUON->RawCluster(ich,1,index1);
2559 if (mRaw1->fClusterType==1 || mRaw1->fClusterType==2) {
2560 Float_t xclust=mCor->fX[3];
2561 Float_t yclust=mCor->fY[3];
2562 AliMUONchamber *iChamber=&(MUON->Chamber(ich));
2563 AliMUONsegmentation *seg = iChamber->GetSegmentationModel(1);
2565 seg->GetPadIxy(xclust,yclust,ix,iy);
2566 Int_t isec=seg->Sector(ix,iy);
2567 printf("nev, CORRELATION with pure background in chamber sector %d %d %d !!!!!!!!!!!!!!!!!!!!!\n",nev,ich+1,isec);
2570 } // end if cluster type on cathode 1
2572 xtrg[maxidg] = (Double_t) mCor->fY[3];
2573 ytrg[maxidg] = (Double_t) mCor->fX[3];
2575 izch[maxidg] = ich+1;
2578 clsize1[maxidg] = mult1;
2579 clsize2[maxidg] = mult2;
2581 cx[maxidg] = 0; // Px/P of hit
2582 cy[maxidg] = 0; // Py/P of hit
2583 cz[maxidg] = 0; // Pz/P of hit
2584 itypg[maxidg] = 0; // particle number
2585 ptotg[maxidg] = 0; // P of hit
2587 pvert1g[maxidg] = 0; // Px vertex
2588 pvert2g[maxidg] = 0; // Py vertex
2589 pvert3g[maxidg] = 0; // Pz vertex
2590 zvertg[maxidg] = 0; // z vertex
2595 } else if (mRaw->fClusterType ==1 && ftrack1 < 0) // background + resonance
2598 // get indexmap and loop over digits to find the signal
2599 Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
2600 gAlice->ResetDigits();
2602 //gAlice->TreeD()->GetEvent(2); // cathode 2
2603 gAlice->TreeD()->GetEvent(nent-1); // cathode 2
2605 //gAlice->TreeD()->GetEvent(1); // cathode 1
2606 gAlice->TreeD()->GetEvent(nent-2); // cathode 1
2609 TClonesArray *MUONdigits = MUON->DigitsAddress(ich);
2610 Int_t mul=mRaw->fMultiplicity;
2612 for (int i=0;i<mul;i++) {
2613 Int_t idx=mRaw->fIndexMap[i];
2614 AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx);
2615 trsign=dig->fTracks[0];
2617 if (trsign > 0 && ihit >= 0) break;
2619 } // loop over indexmap
2621 //printf("trsign, ihit %d %d\n",trsign,ihit);
2622 //printf("signal+background : trsign %d\n",trsign);
2624 if (trsign < 0 || ihit < 0) { // no signal muon was found
2626 if (maxidg<=20000) {
2628 xtrg[maxidg] = (Double_t) mCor->fY[3];
2629 ytrg[maxidg] = (Double_t) mCor->fX[0];
2631 xtrg[maxidg] = (Double_t) mCor->fY[3];
2632 ytrg[maxidg] = (Double_t) mCor->fX[3];
2635 izch[maxidg] = ich+1;
2637 // initialisation of informations which
2638 // can't be reached for background
2640 xgeant[maxidg] = 0; // only for resonances
2641 ygeant[maxidg] = 0; // only for resonances
2642 clsize1[maxidg] = mult1;
2643 clsize2[maxidg] = mult2;
2645 cx[maxidg] = 0; // Px/P of hit
2646 cy[maxidg] = 0; // Py/P of hit
2647 cz[maxidg] = 0; // Pz/P of hit
2648 itypg[maxidg] = 0; // particle number
2649 ptotg[maxidg] = 0; // P of hit
2651 pvert1g[maxidg] = 0; // Px vertex
2652 pvert2g[maxidg] = 0; // Py vertex
2653 pvert3g[maxidg] = 0; // Pz vertex
2658 } else { // signal muon - retrieve info
2659 //printf("inside trsign, ihit %d %d\n",trsign,ihit);
2660 if (maxidg<=20000) {
2662 xtrg[maxidg] = (Double_t) mCor->fY[3];
2663 ytrg[maxidg] = (Double_t) mCor->fX[0];
2665 xtrg[maxidg] = (Double_t) mCor->fY[3];
2666 ytrg[maxidg] = (Double_t) mCor->fX[3];
2668 izch[maxidg] = ich+1;
2669 clsize1[maxidg] = mult1;
2670 clsize2[maxidg] = mult2;
2672 // initialise and set to the correct values
2675 xgeant[maxidg] = 0; // only for resonances
2676 ygeant[maxidg] = 0; // only for resonances
2678 cx[maxidg] = 0; // Px/P of hit
2679 cy[maxidg] = 0; // Py/P of hit
2680 cz[maxidg] = 0; // Pz/P of hit
2681 itypg[maxidg] = 0; // particle number
2682 ptotg[maxidg] = 0; // P of hit
2684 pvert1g[maxidg] = 0; // Px vertex
2685 pvert2g[maxidg] = 0; // Py vertex
2686 pvert3g[maxidg] = 0; // Pz vertex
2688 // try to retrieve info about signal muons
2689 gAlice->ResetHits();
2690 TH->GetEvent(trsign);
2692 TClonesArray *MUONhits = MUON->Hits();
2693 AliMUONhit *mHit= (AliMUONhit*)MUONhits->
2695 TClonesArray *fPartArray = gAlice->Particles();
2697 Int_t nch=mHit->fChamber-1;
2698 //printf("sig+bgr ich, nch %d %d \n",ich,nch);
2700 Int_t ftrack = mHit->fTrack;
2701 Int_t id = ((TParticle*) fPartArray->
2702 UncheckedAt(ftrack))->GetPdgCode();
2703 if (id==kMuonPlus||id==kMuonMinus) {
2704 xgeant[maxidg] = (Double_t) mHit->fY;
2705 ygeant[maxidg] = (Double_t) mHit->fX;
2706 cx[maxidg] = (Double_t) mHit->fCyHit;
2707 cy[maxidg] = (Double_t) mHit->fCxHit;
2708 cz[maxidg] = (Double_t) mHit->fCzHit;
2710 if (id==kMuonPlus) {
2712 } else if (id==kMuonMinus) {
2714 } else itypg[maxidg] = 0;
2716 ptotg[maxidg] = (Double_t) mHit->fPTot;
2717 Part = (TParticle*) fPartArray->
2718 UncheckedAt(ftrack);
2719 Int_t iparent = Part->GetFirstMother();
2721 id2 = ((TParticle*) fPartArray->
2722 UncheckedAt(ftrack))->GetPdgCode();
2727 ip=((TParticle*) fPartArray->
2728 UncheckedAt(iparent))->GetFirstMother();
2730 id2 = ((TParticle*) fPartArray->
2731 UncheckedAt(iparent))->GetPdgCode();
2735 id2 = ((TParticle*) fPartArray->
2736 UncheckedAt(iparent))->GetPdgCode();
2741 Int_t idum = trsign+1;
2743 if (id2==443 || id2==553) {
2745 if (id2==443) id2=114;
2749 idg[maxidg] = 30000*id1+10000*idum+id2;
2750 pvert1g[maxidg] = (Double_t) Part->Py();
2751 pvert2g[maxidg] = (Double_t) Part->Px();
2752 pvert3g[maxidg] = (Double_t) Part->Pz();
2753 zvertg[maxidg] = (Double_t) Part->Vz();
2758 } // sign+bgr, highest bgr
2760 //pure resonance or mixed cluster with the highest
2761 //contribution coming from resonance
2762 if (mRaw->fClusterType >= 1 && ftrack1>=0)
2764 if (maxidg<=20000) {
2766 xtrg[maxidg] = (Double_t) mCor->fY[3];
2767 ytrg[maxidg] = (Double_t) mCor->fX[0];
2769 xtrg[maxidg] = (Double_t) mCor->fY[3];
2770 ytrg[maxidg] = (Double_t) mCor->fX[3];
2772 clsize1[maxidg] = mult1;
2773 clsize2[maxidg] = mult2;
2774 izch[maxidg] = ich+1;
2776 Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
2777 gAlice->ResetDigits();
2779 //gAlice->TreeD()->GetEvent(2); // cathode 2
2780 gAlice->TreeD()->GetEvent(nent-1); // cathode 2
2782 //gAlice->TreeD()->GetEvent(1); // cathode 1
2783 gAlice->TreeD()->GetEvent(nent-2); // cathode 1
2786 TClonesArray *MUONdigits = MUON->DigitsAddress(ich);
2787 Int_t mul=mRaw->fMultiplicity;
2788 for (int i=0;i<mul;i++) {
2789 Int_t idx=mRaw->fIndexMap[i];
2790 AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx);
2792 if (ihit >= 0) break;
2794 } // loop over indexmap
2795 //printf("fClusterType, ihit %d %d \n",mRaw->fClusterType,ihit);
2797 xgeant[maxidg] = 0; // only for resonances
2798 ygeant[maxidg] = 0; // only for resonances
2800 cx[maxidg] = 0; // Px/P of hit
2801 cy[maxidg] = 0; // Py/P of hit
2802 cz[maxidg] = 0; // Pz/P of hit
2803 itypg[maxidg] = 0; // particle number
2804 ptotg[maxidg] = 0; // P of hit
2806 pvert1g[maxidg] = 0; // Px vertex
2807 pvert2g[maxidg] = 0; // Py vertex
2808 pvert3g[maxidg] = 0; // Pz vertex
2811 gAlice->ResetHits();
2812 TH->GetEvent(ftrack1);
2813 TClonesArray *MUONhits = MUON->Hits();
2814 AliMUONhit *mHit= (AliMUONhit*)MUONhits->
2816 TClonesArray *fPartArray = gAlice->Particles();
2818 Int_t nch=mHit->fChamber-1;
2819 //printf("signal ich, nch %d %d \n",ich,nch);
2821 Int_t ftrack = mHit->fTrack;
2822 Int_t id = ((TParticle*) fPartArray->
2823 UncheckedAt(ftrack))->GetPdgCode();
2824 //printf("id %d \n",id);
2825 if (id==kMuonPlus||id==kMuonMinus) {
2826 xgeant[maxidg] = (Double_t) mHit->fY;
2827 ygeant[maxidg] = (Double_t) mHit->fX;
2828 cx[maxidg] = (Double_t) mHit->fCyHit;
2829 cy[maxidg] = (Double_t) mHit->fCxHit;
2830 cz[maxidg] = (Double_t) mHit->fCzHit;
2832 if (id==kMuonPlus) {
2834 } else if (id==kMuonMinus) {
2836 } else itypg[maxidg] = 0;
2838 ptotg[maxidg] = (Double_t) mHit->fPTot;
2839 Part = (TParticle*) fPartArray->
2840 UncheckedAt(ftrack);
2841 Int_t iparent = Part->GetFirstMother();
2843 id2 = ((TParticle*) fPartArray->
2844 UncheckedAt(ftrack))->GetPdgCode();
2849 ip=((TParticle*) fPartArray->
2850 UncheckedAt(iparent))->GetFirstMother();
2852 id2 = ((TParticle*) fPartArray->
2853 UncheckedAt(iparent))->GetPdgCode();
2857 id2 = ((TParticle*) fPartArray->
2858 UncheckedAt(iparent))->GetPdgCode();
2863 Int_t idum = ftrack1+1;
2865 if (id2==443 || id2==553) {
2867 if (id2==443) id2=114;
2870 // printf("id2 %d\n",id2);
2871 idg[maxidg] = 30000*id1+10000*idum+id2;
2872 pvert1g[maxidg] = (Double_t) Part->Py();
2873 pvert2g[maxidg] = (Double_t) Part->Px();
2874 pvert3g[maxidg] = (Double_t) Part->Pz();
2875 zvertg[maxidg] = (Double_t) Part->Vz();
2881 } // if cluster type
2888 cout<<"evenement "<<ievr<<endl;
2890 cout<<"nhittot1="<<nhittot1<<endl;
2892 static Int_t nbres=0;
2893 static Int_t nbcor=0;
2894 if (nres>=19) nbres++;
2895 printf("nres ,nncor - %d %d\n",nres,nncor);
2896 printf("nbres - %d\n",nbres);
2897 if (nncor>=20) nbcor++;
2898 printf("nbcor - %d\n",nbcor);
2899 printf("nuncor - %d\n",nuncor);
2900 printf("nbadcor - %d\n",nbadcor);
2908 void 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)
2911 // Fit a track candidate with the following input parameters:
2912 // INPUT : IVERTEX : vertex flag, if IVERTEX=1 (XVERT,YVERT) are free paramaters
2913 // if IVERTEX=1 (XVERT,YVERT)=(0.,0.)
2914 // PEST(5) : starting value of parameters (minuit)
2915 // PSTEP(5) : step size for parameters (minuit)
2916 // OUTPUT : PXZINV,TPHI,TALAM,XVERT,YVERT : fitted value of the parameters
2918 static Double_t arglist[10];
2919 static Double_t c[5] = {0.4, 0.45, 0.45, 90., 90.};
2920 static Double_t b1, b2, epxz, efi, exs, exvert, eyvert;
2924 TMinuit *gMinuit = new TMinuit(5);
2925 gMinuit->mninit(5,10,7);
2926 gMinuit->SetFCN(fcnf); // constant m.f.
2930 gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
2931 // gMinuit->mnseti('track fitting');
2933 gMinuit->mnparm(0, "invmom", pest[0], pstep[0], -c[0], c[0], ierflg);
2934 gMinuit->mnparm(1, "azimuth", pest[1], pstep[1], -c[1], c[1], ierflg);
2935 gMinuit->mnparm(2, "deep", pest[2], pstep[2], -c[2], c[2], ierflg);
2937 gMinuit->mnparm(3, "x ", pest[3], pstep[3], -c[3], c[3], ierflg);
2938 gMinuit->mnparm(4, "y ", pest[4], pstep[4], -c[4], c[4], ierflg);
2941 gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg);
2942 gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg);
2943 gMinuit->mnexcm("EXIT" , arglist, 0, ierflg);
2945 gMinuit->mnpout(0, chname, pxzinv, epxz, b1, b2, ierflg);
2946 gMinuit->mnpout(1, chname, tphi, efi, b1, b2, ierflg);
2947 gMinuit->mnpout(2, chname, talam, exs, b1, b2, ierflg);
2949 gMinuit->mnpout(3, chname, xvert, exvert, b1, b2, ierflg);
2950 gMinuit->mnpout(4, chname, yvert, eyvert, b1, b2, ierflg);
2957 void fcnf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *pest, Int_t iflag)
2960 // function called by trackf_fit
2962 fcn(npar,grad,fval,pest,iflag,futil);
2965 void 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)
2968 // minuit fits for tracking finding
2970 static Double_t arglist[10];
2971 static Double_t c1[5] = {0.001, 0.001, 0.001, 1., 1.};
2972 static Double_t c2[5] = {0.5, 0.5, 0.5, 120., 120.};
2973 static Double_t emat[9];
2974 static Double_t b1, b2;
2975 Double_t fmin, fedm, errdef;
2976 Int_t npari, nparx, istat;
2981 TMinuit *gMinuit = new TMinuit(5);
2982 gMinuit->mninit(5,10,7);
2983 gMinuit->SetFCN(fcnfitf);
2986 gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
2988 // gMinuit->mnseti('track fitting');
2990 gMinuit->mnparm(0,"invmom", pxzinv, c1[0], -c2[0], c2[0], ierflg); // 0.003, 0.5
2991 gMinuit->mnparm(1,"azimuth ", fis, c1[1], -c2[1], c2[1], ierflg);
2992 gMinuit->mnparm(2,"deep ", alams, c1[2], -c2[2], c2[2], ierflg);
2993 gMinuit->mnparm(3,"xvert", xvert, c1[3], -c2[3], c2[3], ierflg);
2994 gMinuit->mnparm(4,"yvert", yvert, c1[4], -c2[4], c2[4], ierflg);
2996 gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg);
2998 gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg);
2999 gMinuit->mnexcm("EXIT", arglist, 0, ierflg);
3001 gMinuit->mnpout(0, chname, pxzinvf, epxzinv, b1, b2, ierflg);
3002 gMinuit->mnpout(1, chname, fif, efi, b1, b2, ierflg);
3003 gMinuit->mnpout(2, chname, alf, exs, b1, b2, ierflg);
3004 gMinuit->mnpout(3, chname, xvertf, exvert, b1, b2, ierflg);
3005 gMinuit->mnpout(4, chname, yvertf, eyvert, b1, b2, ierflg);
3007 gMinuit->mnemat(emat, 3);
3008 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
3013 void fcnfitf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *xval, Int_t iflag)
3016 // function called by prec_fit
3018 fcnfit(npar,grad,fval,xval,iflag,futil);
3021 ///////////////////// fin modifs perso //////////////////////
3023 ClassImp(AliMUONcluster)
3025 //___________________________________________
3026 AliMUONcluster::AliMUONcluster(Int_t *clhits)
3028 fHitNumber=clhits[0];
3036 ClassImp(AliMUONdigit)
3037 //_____________________________________________________________________________
3038 AliMUONdigit::AliMUONdigit(Int_t *digits)
3041 // Creates a MUON digit object to be updated
3045 fSignal = digits[2];
3046 fPhysics = digits[3];
3050 //_____________________________________________________________________________
3051 AliMUONdigit::AliMUONdigit(Int_t *tracks, Int_t *charges, Int_t *digits)
3054 // Creates a MUON digit object
3058 fSignal = digits[2];
3059 fPhysics = digits[3];
3061 for(Int_t i=0; i<10; i++) {
3062 fTcharges[i] = charges[i];
3063 fTracks[i] = tracks[i];
3067 AliMUONdigit::~AliMUONdigit()
3072 ClassImp(AliMUONlist)
3074 //____________________________________________________________________________
3075 AliMUONlist::AliMUONlist(Int_t ich, Int_t *digits):
3076 AliMUONdigit(digits)
3079 // Creates a MUON digit list object
3083 fTrackList = new TObjArray;
3087 ClassImp(AliMUONhit)
3089 //___________________________________________
3090 AliMUONhit::AliMUONhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
3091 AliHit(shunt, track)
3102 fPHfirst=(Int_t) hits[8];
3103 fPHlast=(Int_t) hits[9];
3111 ClassImp(AliMUONcorrelation)
3112 //___________________________________________
3113 //_____________________________________________________________________________
3114 AliMUONcorrelation::AliMUONcorrelation(Int_t *idx, Float_t *x, Float_t *y)
3117 // Creates a MUON correlation object
3119 for(Int_t i=0; i<4; i++) {
3120 fCorrelIndex[i] = idx[i];
3125 ClassImp(AliMUONRawCluster)
3126 Int_t AliMUONRawCluster::Compare(TObject *obj)
3129 AliMUONRawCluster *raw=(AliMUONRawCluster *)obj;
3130 Float_t r=GetRadius();
3131 Float_t ro=raw->GetRadius();
3133 else if (r<ro) return -1;
3136 AliMUONRawCluster *raw=(AliMUONRawCluster *)obj;
3140 else if (y<yo) return -1;
3145 Int_t AliMUONRawCluster::
3146 BinarySearch(Float_t y, TArrayF coord, Int_t from, Int_t upto)
3148 // Find object using a binary search. Array must first have been sorted.
3149 // Search can be limited by setting upto to desired index.
3151 Int_t low=from, high=upto-1, half;
3154 if(y>coord[half]) low=half;
3160 void AliMUONRawCluster::SortMin(Int_t *idx,Float_t *xdarray,Float_t *xarray,Float_t *yarray,Float_t *qarray, Int_t ntr)
3163 // Get the 3 closest points(cog) one can find on the second cathode
3164 // starting from a given cog on first cathode
3168 // Loop over deltax, only 3 times
3173 Int_t id[3] = {-2,-2,-2};
3174 Float_t jx[3] = {0.,0.,0.};
3175 Float_t jy[3] = {0.,0.,0.};
3176 Float_t jq[3] = {0.,0.,0.};
3177 Int_t jid[3] = {-2,-2,-2};
3180 if (ntr<3) imax=ntr;
3182 for(i=0;i<imax;i++){
3187 if ((i == 1 && j == id[i-1])
3188 ||(i == 2 && (j == id[i-1] || j == id[i-2]))) continue;
3189 if (TMath::Abs(xdarray[j]) < xmin) {
3190 xmin = TMath::Abs(xdarray[j]);
3194 if (xmin != 1001.) {
3221 Int_t AliMUONRawCluster::PhysicsContribution()
3226 for (Int_t i=0; i<fMultiplicity; i++) {
3227 if (fPhysicsMap[i]==2) iPhys++;
3228 if (fPhysicsMap[i]==1) iMixed++;
3229 if (fPhysicsMap[i]==0) iBg++;
3231 if (iMixed==0 && iBg==0) {
3233 } else if ((iPhys != 0 && iBg !=0) || iMixed != 0) {
3241 ClassImp(AliMUONreccluster)
3242 ClassImp(AliMUONsegmentation)
3243 ClassImp(AliMUONresponse)