////////////////////////////////////////////////
#include <TTUBE.h>
+#include <TBRIK.h>
+#include <TRotMatrix.h>
#include <TNode.h>
+#include <TTree.h>
#include <TRandom.h>
#include <TObject.h>
#include <TVector.h>
#include <TObjArray.h>
+#include <TMinuit.h>
+#include <TParticle.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TNtuple.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TDirectory.h>
+#include <TObjectTable.h>
+#include <AliPDG.h>
#include "AliMUON.h"
+#include "TTUBE.h"
+#include "AliMUONClusterFinder.h"
#include "AliRun.h"
#include "AliMC.h"
#include "iostream.h"
#include "AliCallf77.h"
+#ifndef WIN32
+# define reco_init reco_init_
+# define cutpxz cutpxz_
+# define sigmacut sigmacut_
+# define xpreci xpreci_
+# define ypreci ypreci_
+# define reconstmuon reconstmuon_
+# define trackf_read_geant trackf_read_geant_
+# define trackf_read_spoint trackf_read_spoint_
+# define chfill chfill_
+# define chfill2 chfill2_
+# define chf1 chf1_
+# define chfnt chfnt_
+# define hist_create hist_create_
+# define hist_closed hist_closed_
+# define rndm rndm_
+# define fcn fcn_
+# define trackf_fit trackf_fit_
+# define prec_fit prec_fit_
+# define fcnfit fcnfit_
+# define reco_term reco_term_
+#else
+# define reco_init RECO_INIT
+# define cutpxz CUTPXZ
+# define sigmacut SIGMACUT
+# define xpreci XPRECI
+# define ypreci YPRECI
+# define reconstmuon RECONSTMUON
+# define trackf_read_geant TRACKF_READ_GEANT
+# define trackf_read_spoint TRACKF_READ_SPOINT
+# define chfill CHFILL
+# define chfill2 CHFILL2
+# define chf1 CHF1
+# define chfnt CHFNT
+# define hist_create HIST_CREATE
+# define hist_closed HIST_CLOSED
+# define rndm RNDM
+# define fcn FCN
+# define trackf_fit TRACKF_FIT
+# define prec_fit PREC_FIT
+# define fcnfit FCNFIT
+# define reco_term RECO_TERM
+#endif
+
+extern "C"
+{
+void type_of_call reco_init(Double_t &, Double_t &, Double_t &);
+void type_of_call reco_term();
+void type_of_call cutpxz(Double_t &);
+void type_of_call sigmacut(Double_t &);
+void type_of_call xpreci(Double_t &);
+void type_of_call ypreci(Double_t &);
+void type_of_call reconstmuon(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
+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 *);
+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 *);
+void type_of_call chfill(Int_t &, Float_t &, Float_t &, Float_t &);
+void type_of_call chfill2(Int_t &, Float_t &, Float_t &, Float_t &);
+void type_of_call chf1(Int_t &, Float_t &, Float_t &);
+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 *);
+void type_of_call hist_create();
+void type_of_call hist_closed();
+void type_of_call fcnf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
+void type_of_call fcn(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &);
+void type_of_call trackf_fit(Int_t &, Double_t *, Double_t *, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
+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 &);
+void type_of_call fcnfitf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
+void type_of_call fcnfit(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &);
+Float_t type_of_call rndm() {return gRandom->Rndm();}
+}
+
// Static variables for the pad-hit iterator routines
static Int_t sMaxIterPad=0;
static Int_t sCurIterPad=0;
-
+static TTree *TrH1;
+static TTree *TK1;
+static TClonesArray *fHits2; //Listof hits for one track only
+static TClonesArray *fClusters2; //List of clusters for one track only
+static TClonesArray *fParticles2; //List of particles in the Kine tree
ClassImp(AliMUON)
-
//___________________________________________
AliMUON::AliMUON()
{
fClusters = 0;
fNclusters = 0;
fDchambers = 0;
- fRecClusters= 0;
fNdch = 0;
+ fRawClusters= 0;
+ fNrawch = 0;
+ fCathCorrel= 0;
+ fNcorch = 0;
+ fTreeC = 0;
+
+ // modifs perso
+ fSPxzCut = 0;
+ fSSigmaCut = 0;
+ fSXPrec = 0;
+ fSYPrec = 0;
}
//___________________________________________
*/
//End_Html
- fHits = new TClonesArray("AliMUONhit",1000 );
+ fHits = new TClonesArray("AliMUONhit",1000);
fClusters = new TClonesArray("AliMUONcluster",10000);
fNclusters = 0;
fIshunt = 0;
- fNdch = new Int_t[11];
+ fNdch = new Int_t[10];
- fDchambers = new TObjArray(11);
+ fDchambers = new TObjArray(10);
Int_t i;
- for (i=0; i<11 ;i++) {
+ for (i=0; i<10 ;i++) {
(*fDchambers)[i] = new TClonesArray("AliMUONdigit",10000);
fNdch[i]=0;
}
- fRecClusters=new TObjArray(20);
- for (i=0; i<20;i++)
- (*fRecClusters)[i] = new TObjArray(1000);
+ fNrawch = new Int_t[10];
+
+ fRawClusters = new TObjArray(10);
+
+ for (i=0; i<10 ;i++) {
+ (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000);
+ fNrawch[i]=0;
+ }
+
+ fNcorch = new Int_t[10];
+ fCathCorrel = new TObjArray(10);
+ for (i=0; i<10 ;i++) {
+ (*fCathCorrel)[i] = new TClonesArray("AliMUONcorrelation",1000);
+ fNcorch[i]=0;
+ }
+
+ fTreeC = 0;
//
// Transport angular cut
fAccMin=2;
fAccMax=9;
+ // modifs perso
+ fSPxzCut = 3.0;
+ fSSigmaCut = 1.0;
+ fSXPrec = 0.01;
+ fSYPrec = 0.144;
+
SetMarkerColor(kRed);
}
//___________________________________________
AliMUON::~AliMUON()
{
+
+ printf("Calling AliMUON destructor !!!\n");
+
+ Int_t i;
fIshunt = 0;
delete fHits;
delete fClusters;
-// for (Int_t i=0;i<10;i++) {
-// delete (*fDchambers)[i];
-// fNdch[i]=0;
-// }
-// delete fDchambers;
- for (Int_t i=0;i<20;i++)
- delete (*fRecClusters)[i];
- delete fRecClusters;
+ delete fTreeC;
+
+ for (i=0;i<10;i++) {
+ delete (*fDchambers)[i];
+ fNdch[i]=0;
+ }
+ delete fDchambers;
+
+ for (i=0;i<10;i++) {
+ delete (*fRawClusters)[i];
+ fNrawch[i]=0;
+ }
+ delete fRawClusters;
+ for (i=0;i<10;i++) {
+ delete (*fCathCorrel)[i];
+ fNcorch[i]=0;
+ }
+ delete fCathCorrel;
}
//___________________________________________
new(ldigits[fNdch[id]++]) AliMUONdigit(tracks,charges,digits);
}
+//_____________________________________________________________________________
+void AliMUON::AddRawCluster(Int_t id, const AliMUONRawCluster& c)
+{
+ //
+ // Add a MUON digit to the list
+ //
+ TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
+ new(lrawcl[fNrawch[id]++]) AliMUONRawCluster(c);
+}
//_____________________________________________________________________________
-void AliMUON::AddRecCluster(Int_t iCh, Int_t iCat, AliMUONRecCluster* Cluster)
+void AliMUON::AddCathCorrel(Int_t id, Int_t *idx, Float_t *x, Float_t *y)
{
//
- // Add a MUON reconstructed cluster to the list
+ // Add a MUON digit to the list
//
- TObjArray* ClusterList = RecClusters(iCh,iCat);
- ClusterList->Add(Cluster);
+
+ TClonesArray &lcorrel = *((TClonesArray*)(*fCathCorrel)[id]);
+ new(lcorrel[fNcorch[id]++]) AliMUONcorrelation(idx,x,y);
}
//___________________________________________
void AliMUON::BuildGeometry()
{
- TNode *Node, *Top;
- const int kColorMUON = kBlue;
- //
- Top=gAlice->GetGeometry()->GetNode("alice");
-
- // MUON
- const float cz[10] = { 511, 519, 686, 694, 971, 979, 1245, 1253, 1445, 1453};
- const float acc_min = TMath::Tan( 2*.0174532925199432958);
- const float acc_max = TMath::Tan(9*.0174532925199432958);
- float rmin, rmax;
-
- // Chamber 1
- rmin = (cz[0]+0.25)*acc_min;
- rmax = cz[0]*acc_max;
- new TTUBE("S_MUON1","MUON chamber 1","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON1","MUON chamber 1","S_MUON1",0,0,cz[0],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
- // Chamber 2
- rmin = (cz[1]+0.25)*acc_min;
- rmax = cz[1]*acc_max;
- new TTUBE("S_MUON2","MUON chamber 2","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON2","MUON chamber 2","S_MUON2",0,0,cz[1],"");
- fNodes->Add(Node);
- Node->SetLineColor(kColorMUON);
-
- // Chamber 3
- rmin = (cz[2]+0.25)*acc_min;
- rmax = cz[2]*acc_max;
- new TTUBE("S_MUON3","MUON chamber 3","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON3","MUON chamber 3","S_MUON3",0,0,cz[2],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
- // Chamber 4
- rmin = (cz[3]+0.25)*acc_min;
- rmax = cz[3]*acc_max;
- new TTUBE("S_MUON4","MUON chamber 4","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON4","MUON chamber 4","S_MUON4",0,0,cz[3],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
- // Chamber 5
- rmin = 30;
- rmax = cz[4]*acc_max;
- new TTUBE("S_MUON5","MUON chamber 5","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON5","MUON chamber 5","S_MUON5",0,0,cz[4],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
- // Chamber 6
- rmin = 30;
- rmax = cz[5]*acc_max;
- new TTUBE("S_MUON6","MUON chamber 6","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON6","MUON chamber 6","S_MUON6",0,0,cz[5],"");
- fNodes->Add(Node);
- Node->SetLineColor(kColorMUON);
-
- // Chamber 7
- rmin = 30;
- rmax = cz[6]*acc_max;
- new TTUBE("S_MUON7","MUON chamber 7","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON7","MUON chamber 7","S_MUON7",0,0,cz[6],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
- // Chamber 8
- rmin = 30;
- rmax = cz[7]*acc_max;
- new TTUBE("S_MUON8","MUON chamber 8","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON8","MUON chamber 8","S_MUON8",0,0,cz[7],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
- // Chamber 9
- rmin = 30;
- rmax = cz[8]*acc_max;
- new TTUBE("S_MUON9","MUON chamber 9","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON9","MUON chamber 9","S_MUON9",0,0,cz[8],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
- // Chamber 10
- rmin = 30;
- rmax = cz[9]*acc_max;
- new TTUBE("S_MUON10","MUON chamber 10","void",rmin,rmax,0.25);
- Top->cd();
- Node = new TNode("MUON10","MUON chamber 10","S_MUON10",0,0,cz[9],"");
- Node->SetLineColor(kColorMUON);
- fNodes->Add(Node);
-
+ TNode *Node, *NodeF, *Top;
+ const int kColorMUON = kBlue;
+ //
+ Top=gAlice->GetGeometry()->GetNode("alice");
+// MUON
+//
+// z-Positions of Chambers
+ const Float_t cz[5]={511., 686., 971., 1245., 1445.};
+//
+// inner diameter
+ const Float_t dmi[5]={ 35., 47., 67., 86., 100.};
+//
+// outer diameter
+ const Float_t dma[5]={183., 245., 346., 520., 520.};
+
+ TRotMatrix* rot000 = new TRotMatrix("Rot000"," ", 90, 0, 90, 90, 0, 0);
+ TRotMatrix* rot090 = new TRotMatrix("Rot090"," ", 90, 90, 90,180, 0, 0);
+ TRotMatrix* rot180 = new TRotMatrix("Rot180"," ", 90,180, 90,270, 0, 0);
+ TRotMatrix* rot270 = new TRotMatrix("Rot270"," ", 90,270, 90, 0, 0, 0);
+
+
+ float rmin, rmax, dx, dy, dz, dr, zpos;
+ float dzc=4.;
+ char NameChamber[9], NameSense[9], NameFrame[9], NameNode[7];
+ for (Int_t i=0; i<5; i++) {
+ for (Int_t j=0; j<2; j++) {
+ Int_t id=2*i+j+1;
+ if (j==0) {
+ zpos=cz[i]-dzc;
+ } else {
+ zpos=cz[i]+dzc;
+ }
+
+
+ sprintf(NameChamber,"C_MUON%d",id);
+ sprintf(NameSense,"S_MUON%d",id);
+ sprintf(NameFrame,"F_MUON%d",id);
+ rmin = dmi[i]/2.-3;
+ rmax = dma[i]/2.+3;
+ new TTUBE(NameChamber,"Mother","void",rmin,rmax,0.25,1.);
+ rmin = dmi[i]/2.;
+ rmax = dma[i]/2.;
+ new TTUBE(NameSense,"Sens. region","void",rmin,rmax,0.25, 1.);
+ dx=(rmax-rmin)/2;
+ dy=3.;
+ dz=0.25;
+ TBRIK* FMUON = new TBRIK(NameFrame,"Frame","void",dx,dy,dz);
+ Top->cd();
+ sprintf(NameNode,"MUON%d",100+id);
+ Node = new TNode(NameNode,"ChamberNode",NameChamber,0,0,zpos,"");
+ Node->SetLineColor(kColorMUON);
+ fNodes->Add(Node);
+ Node->cd();
+ sprintf(NameNode,"MUON%d",200+id);
+ Node = new TNode(NameNode,"Sens. Region Node",NameSense,0,0,0,"");
+ Node->SetLineColor(kColorMUON);
+ fNodes->Add(Node);
+ Node->cd();
+ dr=dx+rmin;
+ sprintf(NameNode,"MUON%d",300+id);
+ NodeF = new TNode(NameNode,"Frame0",FMUON,dr, 0, 0,rot000,"");
+ NodeF->SetLineColor(kColorMUON);
+ fNodes->Add(NodeF);
+ Node->cd();
+ sprintf(NameNode,"MUON%d",400+id);
+ NodeF = new TNode(NameNode,"Frame1",FMUON,0 ,dr,0,rot090,"");
+ NodeF->SetLineColor(kColorMUON);
+ fNodes->Add(NodeF);
+ Node->cd();
+ sprintf(NameNode,"MUON%d",500+id);
+ NodeF = new TNode(NameNode,"Frame2",FMUON,-dr,0,0,rot180,"");
+ NodeF->SetLineColor(kColorMUON);
+ fNodes->Add(NodeF);
+ Node ->cd();
+ sprintf(NameNode,"MUON%d",600+id);
+ NodeF = new TNode(NameNode,"Frame3",FMUON,0,-dr,0,rot270,"");
+ NodeF->SetLineColor(kColorMUON);
+ fNodes->Add(NodeF);
+ }
+ }
}
+
//___________________________________________
Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
{
// Create Tree branches for the MUON.
const Int_t buffersize = 4000;
- char branchname[20];
+ char branchname[30];
sprintf(branchname,"%sCluster",GetName());
AliDetector::MakeBranch(option);
printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
}
}
-// one branch for rec clusters
- for (i=0; i<20; i++) {
- sprintf(branchname,"%sRecClus%d",GetName(),i+1);
- if (fRecClusters && gAlice->TreeD()) {
- gAlice->TreeR()
- ->Branch(branchname,"TObjArray",
- &((*fRecClusters)[i]), buffersize,0);
- printf("Making Branch %s for clusters in chamber %d\n",
- branchname,i+1);
- }
+
+ //printf("Make Branch - TreeR address %p\n",gAlice->TreeR());
+
+// one branch for raw clusters per chamber
+ for (i=0; i<10 ;i++) {
+ sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
+
+ if (fRawClusters && gAlice->TreeR()) {
+ gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), buffersize);
+ printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
+ }
}
+
}
//___________________________________________
void AliMUON::SetTreeAddress()
{
// Set branch address for the Hits and Digits Tree.
- char branchname[20];
+ char branchname[30];
AliDetector::SetTreeAddress();
TBranch *branch;
TTree *treeH = gAlice->TreeH();
TTree *treeD = gAlice->TreeD();
+ TTree *treeR = gAlice->TreeR();
if (treeH) {
if (fClusters) {
}
}
}
+
+ // printf("SetTreeAddress --- treeR address %p \n",treeR);
+
+ if (treeR) {
+ for (int i=0; i<10; i++) {
+ sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
+ if (fRawClusters) {
+ branch = treeR->GetBranch(branchname);
+ if (branch) branch->SetAddress(&((*fRawClusters)[i]));
+ }
+ }
+ }
+
}
//___________________________________________
void AliMUON::ResetHits()
// Reset number of digits and the digits array for this detector
//
for ( int i=0;i<10;i++ ) {
- if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
+ if ((*fDchambers)[i]) ((TClonesArray*)(*fDchambers)[i])->Clear();
if (fNdch) fNdch[i]=0;
}
}
//____________________________________________
-void AliMUON::ResetRecClusters()
+void AliMUON::ResetRawClusters()
{
//
- // Reset the rec clusters
+ // Reset number of raw clusters and the raw clust array for this detector
//
- for ( int i=0;i<20;i++ ) {
- if ((*fRecClusters)[i]) (*fRecClusters)[i]->Clear();
+ for ( int i=0;i<10;i++ ) {
+ if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
+ if (fNrawch) fNrawch[i]=0;
}
}
+//____________________________________________
+void AliMUON::ResetCorrelation()
+{
+ //
+ // Reset number of correl clusters and the correl clust array for
+ // this detector
+ //
+ for ( int i=0;i<10;i++ ) {
+ if ((*fCathCorrel)[i]) ((TClonesArray*)(*fCathCorrel)[i])->Clear();
+ if (fNcorch) fNcorch[i]=0;
+ }
+}
+
//___________________________________________
void AliMUON::SetPADSIZ(Int_t id, Int_t isec, Float_t p1, Float_t p2)
}
//___________________________________________
-void AliMUON::SetMUCHSP(Int_t id, Float_t p1)
+void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
{
Int_t i=2*(id-1);
- ((AliMUONchamber*) (*fChambers)[i])->SetMUCHSP(p1);
- ((AliMUONchamber*) (*fChambers)[i+1])->SetMUCHSP(p1);
+ ((AliMUONchamber*) (*fChambers)[i])->SetChargeSlope(p1);
+ ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSlope(p1);
}
//___________________________________________
-void AliMUON::SetMUSIGM(Int_t id, Float_t p1, Float_t p2)
+void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
{
Int_t i=2*(id-1);
- ((AliMUONchamber*) (*fChambers)[i])->SetMUSIGM(p1,p2);
- ((AliMUONchamber*) (*fChambers)[i+1])->SetMUSIGM(p1,p2);
+ ((AliMUONchamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
+ ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSpread(p1,p2);
}
//___________________________________________
-void AliMUON::SetRSIGM(Int_t id, Float_t p1)
+void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
{
Int_t i=2*(id-1);
- ((AliMUONchamber*) (*fChambers)[i])->SetRSIGM(p1);
- ((AliMUONchamber*) (*fChambers)[i+1])->SetRSIGM(p1);
+ ((AliMUONchamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
+ ((AliMUONchamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
}
//___________________________________________
-void AliMUON::SetMAXADC(Int_t id, Float_t p1)
+void AliMUON::SetMaxAdc(Int_t id, Float_t p1)
{
Int_t i=2*(id-1);
- ((AliMUONchamber*) (*fChambers)[i])->SetMAXADC(p1);
- ((AliMUONchamber*) (*fChambers)[i+1])->SetMAXADC(p1);
+ ((AliMUONchamber*) (*fChambers)[i])->SetMaxAdc(p1);
+ ((AliMUONchamber*) (*fChambers)[i+1])->SetMaxAdc(p1);
}
//___________________________________________
-void AliMUON::SetSMAXAR(Float_t p1)
+void AliMUON::SetMaxStepGas(Float_t p1)
{
fMaxStepGas=p1;
}
//___________________________________________
-void AliMUON::SetSMAXAL(Float_t p1)
+void AliMUON::SetMaxStepAlu(Float_t p1)
{
fMaxStepAlu=p1;
}
//___________________________________________
-void AliMUON::SetDMAXAR(Float_t p1)
+void AliMUON::SetMaxDestepGas(Float_t p1)
{
fMaxDestepGas=p1;
}
//___________________________________________
-void AliMUON::SetDMAXAL(Float_t p1)
+void AliMUON::SetMaxDestepAlu(Float_t p1)
{
fMaxDestepAlu=p1;
}
//___________________________________________
-void AliMUON::SetMUONACC(Bool_t acc, Float_t angmin, Float_t angmax)
+void AliMUON::SetMuonAcc(Bool_t acc, Float_t angmin, Float_t angmax)
{
fAccCut=acc;
fAccMin=angmin;
((AliMUONchamber*) (*fChambers)[id])->ResponseModel(response);
}
+void AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinder *reconst)
+{
+ ((AliMUONchamber*) (*fChambers)[id])->ReconstructionModel(reconst);
+}
+
void AliMUON::SetNsec(Int_t id, Int_t nsec)
{
((AliMUONchamber*) (*fChambers)[id])->SetNsec(nsec);
void AliMUON::StepManager()
{
printf("Dummy version of muon step -- it should never happen!!\n");
+ /*
const Float_t kRaddeg = 180/TMath::Pi();
AliMC* pMC = AliMC::GetMC();
Int_t nsec, ipart;
Float_t x[4], p[4];
- Float_t pt, th0, th1;
+ Float_t pt, th0, th2;
char proc[5];
if(fAccCut) {
if((nsec=pMC->NSecondaries())>0) {
pMC->ProdProcess(proc);
- if((pMC->TrackPid()==113 || pMC->TrackPid()==114) && !strcmp(proc,"DCAY")) {
+ if((pMC->TrackPid()==443 || pMC->TrackPid()==553) && !strcmp(proc,"DCAY")) {
//
// Check angular acceptance
- //* --- and have muons from resonance decays in the wanted window ---
+ // --- and have muons from resonance decays in the wanted window ---
if(nsec != 2) {
printf(" AliMUON::StepManager: Strange resonance Decay into %d particles\n",nsec);
pMC->StopEvent();
pMC->GetSecondary(0,ipart,x,p);
pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
th0 = TMath::ATan2(pt,p[2])*kRaddeg;
- pMC->GetSecondary(1,ipart,x,p);
+ pMC->GetSecondary(1,ipart,x,p);
pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
- th1 = TMath::ATan2(pt,p[2])*kRaddeg;
+ th2 = TMath::ATan2(pt,p[2])*kRaddeg;
if(!(fAccMin < th0 && th0 < fAccMax) ||
- !(fAccMin < th1 && th1 < fAccMax))
+ !(fAccMin < th2 && th2 < fAccMax))
pMC->StopEvent();
}
}
}
}
+ */
}
-void AliMUON::ReconstructClusters()
+
+void AliMUON::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol)
{
//
-// Initialize the necessary correspondance table
+// Calls the charge disintegration method of the current chamber and adds
+// the simulated cluster to the root treee
//
- static const Int_t kMaxNpadx = 600;
- static const Int_t kMaxNpady = 600;
- Int_t elem[kMaxNpadx*2][kMaxNpady*2];
+ Int_t clhits[7];
+ Float_t newclust[6][500];
+ Int_t nnew;
+
+
//
-// Loop on chambers and on cathode planes
+// Integrated pulse height on chamber
+
+
+ clhits[0]=fNhits+1;
//
- for (Int_t ich=0;ich<10;ich++)
- for (Int_t icat=0;icat<2;icat++) {
- //
- // Get ready the current chamber stuff
- //
- AliMUONchamber* iChamber= &(this->Chamber(ich));
- AliMUONsegmentation* segmentation =
- iChamber->GetSegmentationModel(icat+1);
- if (!segmentation)
- continue;
- TClonesArray *MUONdigits = this->DigitsAddress(ich);
- if (MUONdigits == 0)
- continue;
- cout << "Npx " << segmentation->Npx()
- << " Npy " << segmentation->Npy() << endl;
- //
- // Ready the digits
- //
- gAlice->ResetDigits();
- gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
- Int_t ndigits = MUONdigits->GetEntriesFast();
- if (ndigits == 0)
+//
+ ((AliMUONchamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust);
+// printf("\n Add new clusters %d %f \n", nnew, eloss*1.e9);
+ Int_t ic=0;
+
+//
+// Add new clusters
+ for (Int_t i=0; i<nnew; i++) {
+ if (Int_t(newclust[3][i]) > 0) {
+ ic++;
+// Cathode plane
+ clhits[1] = Int_t(newclust[5][i]);
+// Cluster Charge
+ clhits[2] = Int_t(newclust[0][i]);
+// Pad: ix
+ clhits[3] = Int_t(newclust[1][i]);
+// Pad: iy
+ clhits[4] = Int_t(newclust[2][i]);
+// Pad: charge
+ clhits[5] = Int_t(newclust[3][i]);
+// Pad: chamber sector
+ clhits[6] = Int_t(newclust[4][i]);
+
+ AddCluster(clhits);
+ }
+ }
+// printf("\n %d new clusters added", ic);
+}
+
+void AliMUON::Digitise(Int_t nev,Int_t bgr_ev,Option_t *option, Option_t *,Text_t *filename)
+{
+ // keep galice.root for signal and name differently the file for
+ // background when add! otherwise the track info for signal will be lost !
+
+ static Bool_t first=kTRUE;
+// static TTree *TrH1;
+ static TFile *File;
+ char *Add = strstr(option,"Add");
+ //char *listoftracks = strstr(opt,"listoftracks");
+
+ AliMUONchamber* iChamber;
+ AliMUONsegmentation* segmentation;
+
+
+ Int_t trk[50];
+ Int_t chtrk[50];
+ TObjArray *list=new TObjArray;
+ static TClonesArray *p_adr=0;
+ if(!p_adr) p_adr=new TClonesArray("TVector",1000);
+ Int_t digits[5];
+
+ AliMUON *MUON = (AliMUON *) gAlice->GetModule("MUON");
+ AliMUONHitMap * HitMap[10];
+ for (Int_t i=0; i<10; i++) {HitMap[i]=0;}
+ if (Add ) {
+ if(first) {
+ fFileName=filename;
+ cout<<"filename"<<fFileName<<endl;
+ File=new TFile(fFileName);
+ cout<<"I have opened "<<fFileName<<" file "<<endl;
+ fHits2 = new TClonesArray("AliMUONhit",1000 );
+ fClusters2 = new TClonesArray("AliMUONcluster",10000);
+ }
+ first=kFALSE;
+ File->cd();
+ //File->ls();
+ // Get Hits Tree header from file
+ if(fHits2) fHits2->Clear();
+ if(fClusters2) fClusters2->Clear();
+ if(TrH1) delete TrH1;
+ TrH1=0;
+
+ char treeName[20];
+ sprintf(treeName,"TreeH%d",bgr_ev);
+ TrH1 = (TTree*)gDirectory->Get(treeName);
+ //printf("TrH1 %p of treename %s for event %d \n",TrH1,treeName,bgr_ev);
+
+ if (!TrH1) {
+ printf("ERROR: cannot find Hits Tree for event:%d\n",bgr_ev);
+ }
+ // Set branch addresses
+ TBranch *branch;
+ char branchname[20];
+ sprintf(branchname,"%s",GetName());
+ if (TrH1 && fHits2) {
+ branch = TrH1->GetBranch(branchname);
+ if (branch) branch->SetAddress(&fHits2);
+ }
+ if (TrH1 && fClusters2) {
+ branch = TrH1->GetBranch("MUONCluster");
+ if (branch) branch->SetAddress(&fClusters2);
+ }
+// test
+ //Int_t ntracks1 =(Int_t)TrH1->GetEntries();
+ //printf("background - ntracks1 - %d\n",ntracks1);
+ }
+ //
+ // loop over cathodes
+ //
+ AliMUONHitMap* hm;
+ Int_t countadr=0;
+ for (int icat=0; icat<2; icat++) {
+ Int_t counter=0;
+ for (Int_t i =0; i<10; i++) {
+ iChamber=(AliMUONchamber*) (*fChambers)[i];
+ if (iChamber->Nsec()==1 && icat==1) {
continue;
- printf("Found %d digits for cathode %d in chamber %d \n",
- ndigits,icat,ich+1);
- AliMUONdigit *mdig;
- AliMUONRecCluster *Cluster;
- //
- // Build the correspondance table
- //
- memset(elem,0,sizeof(Int_t)*kMaxNpadx*kMaxNpady*4);
- Int_t digit;
- for (digit=0; digit<ndigits; digit++)
- {
- mdig = (AliMUONdigit*)MUONdigits->UncheckedAt(digit);
- elem[kMaxNpadx+mdig->fPadX][kMaxNpady+mdig->fPadY] = digit+1;
- // because default is 0
+ } else {
+ segmentation=iChamber->GetSegmentationModel(icat+1);
}
- //
- // Declare some useful variables
- //
- Int_t Xlist[10];
- Int_t Ylist[10];
- Int_t Nlist;
- Int_t nclust=0;
- //
- // loop over digits
- //
- for (digit=0;digit<ndigits;digit++) {
- mdig = (AliMUONdigit*)MUONdigits->UncheckedAt(digit);
- //
- // if digit still available, start clustering
- //
- if (elem[kMaxNpadx+mdig->fPadX][kMaxNpady+mdig->fPadY]) {
- Cluster = new AliMUONRecCluster(digit, ich, icat);
- elem[kMaxNpadx+mdig->fPadX][kMaxNpady+mdig->fPadY]=0;
- //
- // loop over the current list of digits
- // and look for neighbours
- //
- for(Int_t clusDigit=Cluster->FirstDigitIndex();
- clusDigit!=Cluster->InvalidDigitIndex();
- clusDigit=Cluster->NextDigitIndex()) {
- AliMUONdigit* pDigit=(AliMUONdigit*)MUONdigits
- ->UncheckedAt(clusDigit);
- segmentation->Neighbours(pDigit->fPadX,pDigit->fPadY,
- &Nlist, Xlist, Ylist);
- for (Int_t Ilist=0;Ilist<Nlist;Ilist++) {
- if (elem[kMaxNpadx+Xlist[Ilist]][kMaxNpady
- +Ylist[Ilist]]) {
- //
- // Add the digit at the end and reset the table
- //
- Cluster->AddDigit(elem[kMaxNpadx+Xlist[Ilist]][kMaxNpady+Ylist[Ilist]]-1);
- elem[kMaxNpadx+Xlist[Ilist]][kMaxNpady
- +Ylist[Ilist]] =0;
- } // if elem
- } // for Ilist
- } // for pDigit
- //
- // Store the cluster (good time to do Cluster polishing)
- //
- segmentation->FitXY(Cluster,MUONdigits);
- nclust ++;
- AddRecCluster(ich,icat,Cluster);
+ HitMap[i] = new AliMUONHitMapA1(segmentation, list);
+ }
+ //printf("Start loop over tracks \n");
+//
+// Loop over tracks
+//
+
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks =(Int_t) TH->GetEntries();
+ //printf("signal - ntracks %d\n",ntracks);
+ Int_t nmuon[10]={0,0,0,0,0,0,0,0,0,0};
+ Float_t xhit[10][2];
+ Float_t yhit[10][2];
+
+ for (Int_t track=0; track<ntracks; track++) {
+ gAlice->ResetHits();
+ TH->GetEvent(track);
+
+//
+// Loop over hits
+ for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
+ mHit;
+ mHit=(AliMUONhit*)MUON->NextHit())
+ {
+ Int_t nch = mHit->fChamber-1; // chamber number
+ if (nch >9) continue;
+ iChamber = &(MUON->Chamber(nch));
+ Int_t rmin = (Int_t)iChamber->RInner();
+ Int_t rmax = (Int_t)iChamber->ROuter();
+ // new 17.07.99
+ if (Add) {
+
+ if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) {
+ xhit[nch][nmuon[nch]]=mHit->fX;
+ yhit[nch][nmuon[nch]]=mHit->fY;
+ nmuon[nch]++;
+ if (nmuon[nch] >2) printf("nmuon %d\n",nmuon[nch]);
+
+ }
+ }
+
+
+
+
+//
+// Loop over pad hits
+ for (AliMUONcluster* mPad=
+ (AliMUONcluster*)MUON->FirstPad(mHit,fClusters);
+ mPad;
+ mPad=(AliMUONcluster*)MUON->NextPad(fClusters))
+ {
+ Int_t cathode = mPad->fCathode; // cathode number
+ Int_t ipx = mPad->fPadX; // pad number on X
+ Int_t ipy = mPad->fPadY; // pad number on Y
+ Int_t iqpad = Int_t(mPad->fQpad*kScale);// charge per pad
+// Int_t iqpad = mPad->fQpad; // charge per pad
+//
+//
+
+ if (cathode != (icat+1)) continue;
+ // fill the info array
+ Float_t thex, they;
+ segmentation=iChamber->GetSegmentationModel(cathode);
+ segmentation->GetPadCxy(ipx,ipy,thex,they);
+ Float_t rpad=TMath::Sqrt(thex*thex+they*they);
+ if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
+
+ new((*p_adr)[countadr++]) TVector(2);
+ TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]);
+ trinfo(0)=(Float_t)track;
+ trinfo(1)=(Float_t)iqpad;
+
+ digits[0]=ipx;
+ digits[1]=ipy;
+ digits[2]=iqpad;
+ digits[3]=iqpad;
+ if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) {
+ digits[4]=mPad->fHitNumber;
+ } else digits[4]=-1;
+
+ AliMUONlist* pdigit;
+ // build the list of fired pads and update the info
+ if (!HitMap[nch]->TestHit(ipx, ipy)) {
+
+ list->AddAtAndExpand(
+ new AliMUONlist(nch,digits),counter);
+
+ HitMap[nch]->SetHit(ipx, ipy, counter);
+ counter++;
+ pdigit=(AliMUONlist*)list->At(list->GetLast());
+ // list of tracks
+ TObjArray *trlist=(TObjArray*)pdigit->TrackList();
+ trlist->Add(&trinfo);
+ } else {
+ pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy);
+ // update charge
+ (*pdigit).fSignal+=iqpad;
+ (*pdigit).fPhysics+=iqpad;
+ // update list of tracks
+ TObjArray* trlist=(TObjArray*)pdigit->TrackList();
+ Int_t last_entry=trlist->GetLast();
+ TVector *ptrk_p=(TVector*)trlist->At(last_entry);
+ TVector &ptrk=*ptrk_p;
+ Int_t last_track=Int_t(ptrk(0));
+ Int_t last_charge=Int_t(ptrk(1));
+ if (last_track==track) {
+ last_charge+=iqpad;
+ trlist->RemoveAt(last_entry);
+ trinfo(0)=last_track;
+ trinfo(1)=last_charge;
+ trlist->AddAt(&trinfo,last_entry);
+ } else {
+ trlist->Add(&trinfo);
+ }
+ // check the track list
+ Int_t nptracks=trlist->GetEntriesFast();
+ if (nptracks > 2) {
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pptrk_p=(TVector*)trlist->At(tr);
+ TVector &pptrk=*pptrk_p;
+ trk[tr]=Int_t(pptrk(0));
+ chtrk[tr]=Int_t(pptrk(1));
+ }
+ } // end if nptracks
+ } // end if pdigit
+ } //end loop over clusters
+ } // hit loop
+ } // track loop
+
+ //Int_t nentr1=list->GetEntriesFast();
+ //printf(" \n counter, nentr1 %d %d\n",counter,nentr1);
+
+ // open the file with background
+
+ if (Add ) {
+ ntracks =(Int_t)TrH1->GetEntries();
+ //printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
+ //printf("background - Start loop over tracks \n");
+//
+// Loop over tracks
+//
+ for (Int_t track=0; track<ntracks; track++) {
+
+ if (fHits2) fHits2->Clear();
+ if (fClusters2) fClusters2->Clear();
+
+ TrH1->GetEvent(track);
+//
+// Loop over hits
+ AliMUONhit* mHit;
+ for(int i=0;i<fHits2->GetEntriesFast();++i)
+ {
+ mHit=(AliMUONhit*) (*fHits2)[i];
+ Int_t nch = mHit->fChamber-1; // chamber number
+ if (nch >9) continue;
+ iChamber = &(MUON->Chamber(nch));
+ Int_t rmin = (Int_t)iChamber->RInner();
+ Int_t rmax = (Int_t)iChamber->ROuter();
+ Float_t xbgr=mHit->fX;
+ Float_t ybgr=mHit->fY;
+ Bool_t cond=kFALSE;
+
+ for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) {
+ Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon])
+ +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]);
+ if (dist<100) cond=kTRUE;
+ }
+ if (!cond) continue;
+
+//
+// Loop over pad hits
+ for (AliMUONcluster* mPad=
+ (AliMUONcluster*)MUON->FirstPad(mHit,fClusters2);
+ mPad;
+ mPad=(AliMUONcluster*)MUON->NextPad(fClusters2))
+ {
+
+ Int_t cathode = mPad->fCathode; // cathode number
+ Int_t ipx = mPad->fPadX; // pad number on X
+ Int_t ipy = mPad->fPadY; // pad number on Y
+ Int_t iqpad = Int_t(mPad->fQpad*kScale);// charge per pad
+// Int_t iqpad = mPad->fQpad; // charge per pad
+
+ if (cathode != (icat+1)) continue;
+ //if (!HitMap[nch]->CheckBoundary()) continue;
+ // fill the info array
+ Float_t thex, they;
+ segmentation=iChamber->GetSegmentationModel(cathode);
+ segmentation->GetPadCxy(ipx,ipy,thex,they);
+ Float_t rpad=TMath::Sqrt(thex*thex+they*they);
+ if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
+
+ new((*p_adr)[countadr++]) TVector(2);
+ TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]);
+ trinfo(0)=-1; // tag background
+ trinfo(1)=-1;
+
+ digits[0]=ipx;
+ digits[1]=ipy;
+ digits[2]=iqpad;
+ digits[3]=0;
+ digits[4]=-1;
+
+ AliMUONlist* pdigit;
+ // build the list of fired pads and update the info
+ if (!HitMap[nch]->TestHit(ipx, ipy)) {
+ list->AddAtAndExpand(new AliMUONlist(nch,digits),counter);
+
+ HitMap[nch]->SetHit(ipx, ipy, counter);
+ counter++;
+
+ pdigit=(AliMUONlist*)list->At(list->GetLast());
+ // list of tracks
+ TObjArray *trlist=(TObjArray*)pdigit->
+ TrackList();
+ trlist->Add(&trinfo);
+ } else {
+ pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy);
+ // update charge
+ (*pdigit).fSignal+=iqpad;
+
+ // update list of tracks
+ TObjArray* trlist=(TObjArray*)pdigit->
+ TrackList();
+ Int_t last_entry=trlist->GetLast();
+ TVector *ptrk_p=(TVector*)trlist->
+ At(last_entry);
+ TVector &ptrk=*ptrk_p;
+ Int_t last_track=Int_t(ptrk(0));
+ if (last_track==-1) {
+ continue;
+ } else {
+ trlist->Add(&trinfo);
+ }
+ // check the track list
+ Int_t nptracks=trlist->GetEntriesFast();
+ if (nptracks > 0) {
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pptrk_p=(TVector*)trlist->At(tr);
+ TVector &pptrk=*pptrk_p;
+ trk[tr]=Int_t(pptrk(0));
+ chtrk[tr]=Int_t(pptrk(1));
+ }
+ } // end if nptracks
+ } // end if pdigit
+ } //end loop over clusters
+ } // hit loop
+ } // track loop
+ //Int_t nentr2=list->GetEntriesFast();
+ //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
+ TTree *fAli=gAlice->TreeK();
+ TFile *file;
+
+ if (fAli) file =fAli->GetCurrentFile();
+ file->cd();
+ } // if Add
+
+ Int_t tracks[10];
+ Int_t charges[10];
+ //cout<<"start filling digits \n "<<endl;
+ // const Float_t zero_supm = 6.;
+ Int_t nentries=list->GetEntriesFast();
+ //printf(" \n \n nentries %d \n",nentries);
+ // start filling the digits
+
+ for (Int_t nent=0;nent<nentries;nent++) {
+ AliMUONlist *address=(AliMUONlist*)list->At(nent);
+ if (address==0) continue;
+ Int_t ich=address->fChamber;
+ Int_t q=address->fSignal;
+ iChamber=(AliMUONchamber*) (*fChambers)[ich];
+ AliMUONresponse * response=iChamber->GetResponseModel();
+ Int_t adcmax= (Int_t) response->MaxAdc();
+ // add white noise and do zero-suppression and signal truncation
+ Float_t MeanNoise = gRandom->Gaus(1, 0.2);
+ Float_t Noise = gRandom->Gaus(0, MeanNoise);
+ q+=(Int_t)Noise;
+ if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise;
+ if ( q <= zero_supm ) continue;
+ if ( q > adcmax) q=adcmax;
+ digits[0]=address->fPadX;
+ digits[1]=address->fPadY;
+ digits[2]=q;
+ digits[3]=address->fPhysics;
+ digits[4]=address->fHit;
+ //printf("fSignal, fPhysics fTrack %d %d %d \n",digits[2],digits[3],digits[4]);
+
+ TObjArray* trlist=(TObjArray*)address->TrackList();
+ Int_t nptracks=trlist->GetEntriesFast();
+ //printf("nptracks, trlist %d %p\n",nptracks,trlist);
+
+ // this was changed to accomodate the real number of tracks
+ if (nptracks > 10) {
+ cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
+ nptracks=10;
+ }
+ if (nptracks > 2) {
+ printf("Attention - nptracks > 2 %d \n",nptracks);
+ printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
}
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pp_p=(TVector*)trlist->At(tr);
+ if(!pp_p ) printf("pp_p - %p\n",pp_p);
+ TVector &pp =*pp_p;
+ tracks[tr]=Int_t(pp(0));
+ charges[tr]=Int_t(pp(1));
+ //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]);
+ } //end loop over list of tracks for one pad
+ // Sort list of tracks according to charge
+ if (nptracks > 1) {
+ SortTracks(tracks,charges,nptracks);
+ }
+ if (nptracks < 10 ) {
+ for (Int_t i=nptracks; i<10; i++) {
+ tracks[i]=0;
+ charges[i]=0;
+ }
+ }
+
+ // fill digits
+ MUON->AddDigits(ich,tracks,charges,digits);
+ }
+ //cout<<"I'm out of the loops for digitisation"<<endl;
+ gAlice->TreeD()->Fill();
+ TTree *TD=gAlice->TreeD();
+
+ Stat_t ndig=TD->GetEntries();
+ cout<<"number of digits "<<ndig<<endl;
+ TClonesArray *fDch;
+ for (int k=0;k<10;k++) {
+ fDch= MUON->DigitsAddress(k);
+ int ndig=fDch->GetEntriesFast();
+ printf (" i, ndig %d %d \n",k,ndig);
+ }
+
+ MUON->ResetDigits();
+ list->Delete();
+ for(Int_t ii=0;ii<10;++ii) {
+ if (HitMap[ii]) {
+ hm=HitMap[ii];
+ delete hm;
+ HitMap[ii]=0;
}
- printf("===> %d Clusters\n",nclust);
- } // for icat
+ }
+
+ } //end loop over cathodes
+
+ char hname[30];
+ sprintf(hname,"TreeD%d",nev);
+ gAlice->TreeD()->Write(hname);
+ // reset tree
+ gAlice->TreeD()->Reset();
+ delete list;
+ //Int_t nadr=p_adr->GetEntriesFast();
+ // printf(" \n \n nadr %d \n",nadr);
+
+ p_adr->Clear();
+ // gObjectTable->Print();
+
+}
+
+void AliMUON::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
+{
+ //
+ // Sort the list of tracks contributing to a given digit
+ // Only the 3 most significant tracks are acctually sorted
+ //
+
+ //
+ // Loop over signals, only 3 times
+ //
+
+ Int_t qmax;
+ Int_t jmax;
+ Int_t idx[3] = {-2,-2,-2};
+ Int_t jch[3] = {-2,-2,-2};
+ Int_t jtr[3] = {-2,-2,-2};
+ Int_t i,j,imax;
+
+ if (ntr<3) imax=ntr;
+ else imax=3;
+ for(i=0;i<imax;i++){
+ qmax=0;
+ jmax=0;
+
+ for(j=0;j<ntr;j++){
+
+ if((i == 1 && j == idx[i-1])
+ ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
+
+ if(charges[j] > qmax) {
+ qmax = charges[j];
+ jmax=j;
+ }
+ }
+
+ if(qmax > 0) {
+ idx[i]=jmax;
+ jch[i]=charges[jmax];
+ jtr[i]=tracks[jmax];
+ }
+
+ }
+
+ for(i=0;i<3;i++){
+ if (jtr[i] == -2) {
+ charges[i]=0;
+ tracks[i]=0;
+ } else {
+ charges[i]=jch[i];
+ tracks[i]=jtr[i];
+ }
+ }
+
}
+void AliMUON::FindClusters(Int_t nev,Int_t last_entry)
+{
+
+//
+// Loop on chambers and on cathode planes
+//
+ for (Int_t icat=0;icat<2;icat++) {
+ gAlice->ResetDigits();
+ gAlice->TreeD()->GetEvent(last_entry+icat); // spurious +1 ...
+ if (nev < 10) printf("last_entry , icat - %d %d \n",last_entry,icat);
+ //gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
+
+ for (Int_t ich=0;ich<10;ich++) {
+ AliMUONchamber* iChamber=(AliMUONchamber*) (*fChambers)[ich];
+ TClonesArray *MUONdigits = this->DigitsAddress(ich);
+ if (MUONdigits == 0) continue;
+ //
+ // Get ready the current chamber stuff
+ //
+ AliMUONresponse* response = iChamber->GetResponseModel();
+ AliMUONsegmentation* seg = iChamber->GetSegmentationModel(icat+1);
+ AliMUONClusterFinder* rec = iChamber->GetReconstructionModel();
+ //printf("icat, ich, seg - %d %d %p\n",icat,ich,seg);
+ if (seg) {
+ rec->SetSegmentation(seg);
+ rec->SetResponse(response);
+ rec->SetDigits(MUONdigits);
+ rec->SetChamber(ich);
+ if (nev==0) rec->CalibrateCOG();
+ rec->FindRawClusters();
+ }
+ //printf("Finish FindRawClusters for cathode %d in chamber %d\n",icat,ich);
+
+ TClonesArray *fRch;
+ fRch=RawClustAddress(ich);
+ fRch->Sort();
+ // it seems to work
+
+
+ } // for ich
+ // fill the tree
+ TTree *TR=gAlice->TreeR();
+
+ gAlice->TreeR()->Fill();
+
+ Stat_t nent=TR->GetEntries();
+ cout<<"number of entries "<<nent<<endl;
+ TClonesArray *fRch;
+ for (int i=0;i<10;i++) {
+ fRch=RawClustAddress(i);
+ int nraw=fRch->GetEntriesFast();
+ printf (" i, nraw %d %d \n",i,nraw);
+ }
+ ResetRawClusters();
+
+ } // for icat
+
+ char hname[30];
+ sprintf(hname,"TreeR%d",nev);
+ gAlice->TreeR()->Write(hname);
+ gAlice->TreeR()->Reset();
+
+ //gObjectTable->Print();
+
+}
//______________________________________________________________________________
+//_____________________________________________________________________________
+void AliMUON::CathodeCorrelation(Int_t nev)
+{
+
+// Correlates the clusters on the two cathode planes and build a list of
+// other possible combinations (potential ghosts) - for the moment use the
+// criteria of minimum distance between the CoGs of the two correlated
+// clusters
+
+
+//
+// Loop on chambers and on clusters on the cathode plane with the highest
+// number of clusters
+
+ static Bool_t first=kTRUE;
+
+ AliMUONRawCluster *mRaw1;
+ AliMUONRawCluster *mRaw2;
+ AliMUONchamber *iChamber;
+ AliMUONsegmentation *seg;
+ TArrayF x1, y1, x2, y2, q1, q2;
+ x1.Set(5000);
+ x2.Set(5000);
+ y1.Set(5000);
+ y2.Set(5000);
+ q1.Set(5000);
+ q2.Set(5000);
+
+// Get pointers to Alice detectors and Digits containers
+ TTree *TR = gAlice->TreeR();
+ Int_t nent=(Int_t)TR->GetEntries();
+ if (nev < 10) printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
+
+
+ Int_t idx[4];
+ Float_t xc2[4],yc2[4];
+ Float_t xrec2, yrec2;
+ Float_t xd0, xdif, ydif;
+ Float_t ysrch,xd,xmax,ymax;
+ Int_t ilow, iup, iraw1, i;
+ //
+ Float_t xarray[50];
+ Float_t xdarray[50];
+ Float_t yarray[50];
+ Float_t qarray[50];
+ Int_t idx2[50];
+
+ // Int_t nraw[2], entry,cathode;
+
+ for (i=0;i<50;i++) {
+ xdarray[i]=1100.;
+ xarray[i]=0.;
+ yarray[i]=0.;
+ qarray[i]=0.;
+ idx2[i]=-1;
+ }
+ for (i=0;i<4;i++) {
+ idx[i]=-1;
+ xc2[i]=0.;
+ yc2[i]=0.;
+ }
+
+ // access to the Raw Clusters tree
+ for (Int_t ich=0;ich<10;ich++) {
+ iChamber = &(Chamber(ich));
+ TClonesArray *MUONrawclust = RawClustAddress(ich);
+ ResetRawClusters();
+ TR->GetEvent(nent-2);
+ //TR->GetEvent(1);
+ Int_t nrawcl1 = MUONrawclust->GetEntries();
+ // printf("Found %d raw clusters for cathode 1 in chamber %d \n"
+ // ,nrawcl1,ich+1);
+ if (!nrawcl1) continue;
+
+ seg = iChamber->GetSegmentationModel(1);
+ // loop over raw clusters of first cathode
+ for (iraw1=0; iraw1<nrawcl1; iraw1++) {
+ mRaw1= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw1);
+ x1[iraw1]=mRaw1->fX;
+ y1[iraw1]=mRaw1->fY;
+ q1[iraw1]=(Float_t)mRaw1->fQ; //maybe better fPeakSignal
+ } // rawclusters cathode 1
+//
+ // Get information from 2nd cathode
+ ResetRawClusters();
+ TR->GetEvent(nent-1);
+ //TR->GetEvent(2);
+ Int_t nrawcl2 = MUONrawclust->GetEntries();
+ if (!nrawcl2) {
+ for (iraw1=0; iraw1<nrawcl1; iraw1++) {
+ idx[3]=iraw1;
+ xc2[3]=x1[iraw1];
+ yc2[3]=y1[iraw1];
+ //printf("nrawcl2 is zero - idx[0] %d\n",idx[0]);
+
+ AddCathCorrel(ich,idx,xc2,yc2);
+ // reset
+ idx[3]=-1;
+ xc2[3]=0.;
+ yc2[3]=0.;
+
+ } // store information from cathode 1 only
+ } else {
+ // printf("Found %d raw clusters for cathode 2 in chamber %d \n",
+ // nrawcl2, ich+1);
+
+ for (Int_t iraw2=0; iraw2<nrawcl2; iraw2++) {
+ mRaw2= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw2);
+ x2[iraw2]=mRaw2->fX;
+ y2[iraw2]=mRaw2->fY;
+ q2[iraw2]=(Float_t)mRaw2->fQ;
+ } // rawclusters cathode 2
+//
+// Initalisation finished
+ for (iraw1=0; iraw1<nrawcl1; iraw1++) {
+ // find the sector
+ Int_t ix,iy;
+ seg->GetPadIxy(x1[iraw1],y1[iraw1],ix,iy);
+ Int_t isec=seg->Sector(ix,iy);
+ // range to look for ghosts ?!
+ if (ich < 5) {
+ ymax = seg->Dpy(isec)*7/2;
+ xmax = seg->Dpx(isec)*7/2;
+ } else {
+ ymax = seg->Dpy(isec)*13/2;
+ xmax = seg->Dpx(isec)*3/2;
+ }
+ ysrch=ymax+y1[iraw1];
+
+ ilow = AliMUONRawCluster::
+ BinarySearch(ysrch-2*ymax,y2,0,nrawcl2+1);
+ iup= AliMUONRawCluster::
+ BinarySearch(ysrch,y2,ilow,nrawcl2+1);
+ if (ilow<0 || iup <0 || iup>nrawcl2) continue;
+ Int_t counter=0;
+ for (Int_t iraw2=ilow; iraw2<=iup; iraw2++) {
+ xrec2=x2[iraw2];
+ yrec2=y2[iraw2];
+ xdif=x1[iraw1]-xrec2;
+ ydif=y1[iraw1]-yrec2;
+ xd=TMath::Sqrt(xdif*xdif+ydif*ydif);
+ if (iraw2==ilow) {
+ if (ilow==iup)
+ xd0=TMath::
+ Sqrt(2*xmax*2*xmax+2*ymax*2*ymax);
+ else xd0=101.;
+ }
+ Float_t qdif=TMath::Abs(q1[iraw1]-q2[iraw2])/q1[iraw1];
+
+ if (x1[iraw1]*xrec2 > 0) {
+ if (xd <= xd0 ) {
+// printf("q1, q2 qdif % f %f %f \n",q1[iraw1],q2[iraw2],qdif);
+// printf("x1, x2 y1 y2 % f %f %f %f \n",x1[iraw1],xrec2,y1[iraw1],yrec2);
+ //if (qdif <0.3) { //check this number
+
+ xd0=xd;
+ idx2[counter]=iraw2;
+ xdarray[counter]=xd;
+ xarray[counter]=xdif;
+ yarray[counter]=ydif;
+ qarray[counter]=qdif;
+ counter++;
+ // }
+
+ }
+ } // check for same quadrant
+ } // loop over 2nd cathode range
+
+
+ if (counter >=2) {
+ AliMUONRawCluster::
+ SortMin(idx2,xdarray,xarray,yarray,qarray,counter);
+ if (xdarray[0]<seg->Dpx(isec) && xdarray[1]<seg->Dpx(isec)) {
+ if (qarray[0]>qarray[1]){
+ Int_t swap=idx2[0];
+ idx2[0]=idx2[1];
+ idx2[1]=swap;
+ }
+ }
+ }
+ int imax;
+ if (counter <3) imax=counter;
+ else imax=3;
+
+ for (int i=0;i<imax;i++) {
+ if (idx2[i] >= 0 && idx2[i] < nrawcl2) {
+ if (xarray[i] > xmax || yarray[i] > 2*ymax)
+ continue;
+ idx[i]=idx2[i];
+ xc2[i]=x2[idx2[i]];
+ yc2[i]=y2[idx2[i]];
+ }
+ }
+ // add info about the cluster on the 'starting' cathode
+
+ idx[3]=iraw1;
+ xc2[3]=x1[iraw1];
+ yc2[3]=y1[iraw1];
+ //if (idx[0] <0) printf("iraw1 imax idx2[0] idx[0] %d %d %d %d\n",iraw1,imax,idx2[0],idx[0]);
+ AddCathCorrel(ich,idx,xc2,yc2);
+ // reset
+ for (Int_t ii=0;ii<counter;ii++) {
+ xdarray[ii]=1100.;
+ xarray[ii]=0.;
+ yarray[ii]=0.;
+ qarray[ii]=0.;
+ idx2[ii]=-1;
+ }
+ for (Int_t iii=0;iii<3;iii++) {
+ idx[iii]=-1;
+ xc2[iii]=0.;
+ yc2[iii]=0.;
+ }
+ } // iraw1
+ }
+ x1.Reset();
+ x2.Reset();
+ y1.Reset();
+ y2.Reset();
+ q1.Reset();
+ q2.Reset();
+ } //ich
+//
+ if (first) {
+ MakeTreeC("C");
+ first=kFALSE;
+ }
+ TTree *TC=TreeC();
+ TC->Fill();
+ //Int_t nentries=(Int_t)TC->GetEntries();
+ //cout<<"number entries in tree of correlated clusters "<<nentries<<endl;
+ TClonesArray *fCch;
+ static Int_t countev=0;
+ Int_t countch=0;
+
+ for (Int_t ii=0;ii<10;ii++) {
+ fCch= CathCorrelAddress(ii);
+ Int_t ncor=fCch->GetEntriesFast();
+ printf (" ii, ncor %d %d \n",ii,ncor);
+ if (ncor>=2) countch++;
+ }
+
+ // write
+ char hname[30];
+ sprintf(hname,"TreeC%d",nev);
+ TC->Write(hname);
+ // reset tree
+ ResetCorrelation();
+ TC->Reset();
+
+ if (countch==10) countev++;
+ printf("countev - %d\n",countev);
+
+// gObjectTable->Print();
+
+
+}
+
+
+//_____________________________________________________________________________
+
+void AliMUON::MakeTreeC(Option_t *option)
+{
+ char *C = strstr(option,"C");
+ if (C && !fTreeC) fTreeC = new TTree("TC","CathodeCorrelation");
+
+// Create a branch for correlation
+
+ const Int_t buffersize = 4000;
+ char branchname[30];
+
+// one branch for correlation per chamber
+ for (int i=0; i<10 ;i++) {
+ sprintf(branchname,"%sCorrelation%d",GetName(),i+1);
+
+ if (fCathCorrel && fTreeC) {
+ TreeC()->Branch(branchname,&((*fCathCorrel)[i]), buffersize);
+ printf("Making Branch %s for correlation in chamber %d\n",branchname,i+1);
+ }
+ }
+}
+
+//_____________________________________________________________________________
+void AliMUON::GetTreeC(Int_t event)
+{
+
+ // set the branch address
+ char treeName[20];
+ char branchname[30];
+
+ ResetCorrelation();
+ if (fTreeC) {
+ delete fTreeC;
+ }
+
+ sprintf(treeName,"TreeC%d",event);
+ fTreeC = (TTree*)gDirectory->Get(treeName);
+
+
+ TBranch *branch;
+ if (fTreeC) {
+ for (int i=0; i<10; i++) {
+ sprintf(branchname,"%sCorrelation%d",GetName(),i+1);
+ if (fCathCorrel) {
+ branch = fTreeC->GetBranch(branchname);
+ if (branch) branch->SetAddress(&((*fCathCorrel)[i]));
+ }
+ }
+ } else {
+ printf("ERROR: cannot find CathodeCorrelation Tree for event:%d\n",event);
+ }
+
+ // gObjectTable->Print();
+
+}
+
+
void AliMUON::Streamer(TBuffer &R__b)
{
// Stream an object of class AliMUON.
AliMUONsegmentation *segmentation;
AliMUONresponse *response;
TClonesArray *digitsaddress;
+ TClonesArray *rawcladdress;
+ TClonesArray *corcladdress;
+ // TObjArray *clustaddress;
if (R__b.IsReading()) {
Version_t R__v = R__b.ReadVersion(); if (R__v) { }
R__b >> fNclusters;
R__b >> fClusters; // diff
R__b >> fDchambers;
+ R__b >> fRawClusters;
+ R__b >> fCathCorrel;
R__b.ReadArray(fNdch);
+ R__b.ReadArray(fNrawch);
+ R__b.ReadArray(fNcorch);
//
R__b >> fAccCut;
R__b >> fAccMin;
R__b >> fAccMax;
//
+ // modifs perso
+ R__b >> fSPxzCut;
+ R__b >> fSSigmaCut;
+ R__b >> fSXPrec;
+ R__b >> fSYPrec;
+ //
R__b >> fChambers;
// Stream chamber related information
for (Int_t i =0; i<10; i++) {
response->Streamer(R__b);
digitsaddress=(TClonesArray*) (*fDchambers)[i];
digitsaddress->Streamer(R__b);
+ rawcladdress=(TClonesArray*) (*fRawClusters)[i];
+ rawcladdress->Streamer(R__b);
+ corcladdress=(TClonesArray*) (*fCathCorrel)[i];
+ corcladdress->Streamer(R__b);
}
} else {
R__b << fNclusters;
R__b << fClusters; // diff
R__b << fDchambers;
+ R__b << fRawClusters;
+ R__b << fCathCorrel;
R__b.WriteArray(fNdch, 10);
+ R__b.WriteArray(fNrawch, 10);
+ R__b.WriteArray(fNcorch, 10);
//
R__b << fAccCut;
R__b << fAccMin;
R__b << fAccMax;
//
+ // modifs perso
+ R__b << fSPxzCut;
+ R__b << fSSigmaCut;
+ R__b << fSXPrec;
+ R__b << fSYPrec;
+ //
R__b << fChambers;
// Stream chamber related information
for (Int_t i =0; i<10; i++) {
}
response=iChamber->GetResponseModel();
response->Streamer(R__b);
-
digitsaddress=(TClonesArray*) (*fDchambers)[i];
digitsaddress->Streamer(R__b);
+ rawcladdress=(TClonesArray*) (*fRawClusters)[i];
+ rawcladdress->Streamer(R__b);
+ corcladdress=(TClonesArray*) (*fCathCorrel)[i];
+ corcladdress->Streamer(R__b);
}
}
}
-AliMUONcluster* AliMUON::FirstPad(AliMUONhit* hit)
+AliMUONcluster* AliMUON::FirstPad(AliMUONhit* hit, TClonesArray *clusters)
{
//
// Initialise the pad iterator
// Return the address of the first padhit for hit
- TClonesArray *theClusters = Clusters();
+ TClonesArray *theClusters = clusters;
Int_t nclust = theClusters->GetEntriesFast();
if (nclust && hit->fPHlast > 0) {
sMaxIterPad=hit->fPHlast;
sCurIterPad=hit->fPHfirst;
- return (AliMUONcluster*) fClusters->UncheckedAt(sCurIterPad-1);
+ return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1);
} else {
return 0;
}
}
-AliMUONcluster* AliMUON::NextPad()
+AliMUONcluster* AliMUON::NextPad(TClonesArray *clusters)
{
sCurIterPad++;
if (sCurIterPad <= sMaxIterPad) {
- return (AliMUONcluster*) fClusters->UncheckedAt(sCurIterPad-1);
+ return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1);
} else {
return 0;
}
}
+//////////////////////////// modifs perso ///////////////
+
+static TTree *ntuple_global;
+static TFile *hfile_global;
+
+// variables of the tracking ntuple
+struct {
+ Int_t ievr; // number of event
+ Int_t ntrackr; // number of tracks per event
+ Int_t istatr[500]; // 1 = good muon, 2 = ghost, 0 = something else
+ Int_t isignr[500]; // sign of the track
+ Float_t pxr[500]; // x momentum of the reconstructed track
+ Float_t pyr[500]; // y momentum of the reconstructed track
+ Float_t pzr[500]; // z momentum of the reconstructed track
+ Float_t zvr[500]; // z vertex
+ Float_t chi2r[500]; // chi2 of the fit of the track with the field map
+ Float_t pxv[500]; // x momentum at vertex
+ Float_t pyv[500]; // y momentum at vertex
+ Float_t pzv[500]; // z momentum at vertex
+} ntuple_st;
+
+AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
+{
+ TClonesArray *MUONrawclust = RawClustAddress(ichamber);
+ ResetRawClusters();
+ TTree *TR = gAlice->TreeR();
+ Int_t nent=(Int_t)TR->GetEntries();
+ TR->GetEvent(nent-2+icathod-1);
+ //TR->GetEvent(icathod);
+ //Int_t nrawcl = (Int_t)MUONrawclust->GetEntriesFast();
+
+ AliMUONRawCluster * mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(icluster);
+ //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
+
+ return mRaw;
+}
+
+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)
+{
+ //
+ // open kine and hits tree of background file for reconstruction of geant hits
+ // call tracking fortran program
+ static Bool_t first=kTRUE;
+ static TFile *File;
+ char *Add = strstr(option,"Add");
+
+ if (Add ) { // only in case of background with geant hits
+ if(first) {
+ fFileName=filename;
+ cout<<"filename "<<fFileName<<endl;
+ File=new TFile(fFileName);
+ cout<<"I have opened "<<fFileName<<" file "<<endl;
+ fHits2 = new TClonesArray("AliMUONhit",1000 );
+ fParticles2 = new TClonesArray("GParticle",1000);
+ first=kFALSE;
+ }
+ File->cd();
+ if(fHits2) fHits2->Clear();
+ if(fParticles2) fParticles2->Clear();
+ if(TrH1) delete TrH1;
+ TrH1=0;
+ if(TK1) delete TK1;
+ TK1=0;
+ // Get Hits Tree header from file
+ char treeName[20];
+ sprintf(treeName,"TreeH%d",bgd_ev);
+ TrH1 = (TTree*)gDirectory->Get(treeName);
+ if (!TrH1) {
+ printf("ERROR: cannot find Hits Tree for event:%d\n",bgd_ev);
+ }
+ // set branch addresses
+ TBranch *branch;
+ char branchname[30];
+ sprintf(branchname,"%s",GetName());
+ if (TrH1 && fHits2) {
+ branch = TrH1->GetBranch(branchname);
+ if (branch) branch->SetAddress(&fHits2);
+ }
+ TrH1->GetEntries();
+ // get the Kine tree
+ sprintf(treeName,"TreeK%d",bgd_ev);
+ TK1 = (TTree*)gDirectory->Get(treeName);
+ if (!TK1) {
+ printf("ERROR: cannot find Kine Tree for event:%d\n",bgd_ev);
+ }
+ // set branch addresses
+ if (TK1)
+ TK1->SetBranchAddress("Particles", &fParticles2);
+ TK1->GetEvent(0);
+
+ // get back to the first file
+ TTree *TK = gAlice->TreeK();
+ TFile *file1 = 0;
+ if (TK) file1 = TK->GetCurrentFile();
+ file1->cd();
+
+ } // end if Add
+
+ // call tracking fortran program
+ reconstmuon(ifit,idebug,nev,idres,ireadgeant);
+}
+
+
+void AliMUON::InitTracking(Double_t &seff, Double_t &sb0, Double_t &sbl3)
+{
+ //
+ // introduce in fortran program somme parameters and cuts for tracking
+ // create output file "reconst.root" (histos + ntuple)
+ cutpxz(fSPxzCut); // Pxz cut (GeV/c) to begin the track finding
+ sigmacut(fSSigmaCut); // Number of sigmas delimiting the searching areas
+ xpreci(fSXPrec); // Chamber precision in X (cm)
+ ypreci(fSYPrec); // Chamber precision in Y (cm)
+ reco_init(seff,sb0,sbl3);
+}
+
+void AliMUON::FinishEvent()
+{
+ TTree *TK = gAlice->TreeK();
+ TFile *file1 = 0;
+ if (TK) file1 = TK->GetCurrentFile();
+ file1->cd();
+}
+
+void AliMUON::CloseTracking()
+{
+ //
+ // write histos and ntuple to "reconst.root" file
+ reco_term();
+}
+
+void chfill(Int_t &id, Float_t &x, Float_t &, Float_t &)
+{
+ //
+ // fill histo like hfill in fortran
+ char name[5];
+ sprintf(name,"h%d",id);
+ TH1F *h1 = (TH1F*) gDirectory->Get(name);
+ h1->Fill(x);
+}
+
+void chfill2(Int_t &id, Float_t &x, Float_t &y, Float_t &w)
+{
+ //
+ // fill histo like hfill2 in fortran
+ char name[5];
+ sprintf(name,"h%d",id);
+ TH2F *h2 = (TH2F*) gDirectory->Get(name);
+ h2->Fill(x,y,w);
+}
+
+void chf1(Int_t &id, Float_t &x, Float_t &w)
+{
+ //
+ // fill histo like hf1 in fortran
+ char name[5];
+ sprintf(name,"h%d",id);
+ TH1F *h1 = (TH1F*) gDirectory->Get(name);
+ h1->Fill(x,w);
+}
+
+void hist_create()
+{
+ //
+ // Create an output file ("reconst.root")
+ // Create some histograms and an ntuple
+
+ hfile_global = new TFile("reconst.root","RECREATE","Ntuple - reconstruction");
+
+ ntuple_global = new TTree("ntuple","Reconst ntuple");
+ ntuple_global->Branch("ievr",&ntuple_st.ievr,"ievr/I");
+ ntuple_global->Branch("ntrackr",&ntuple_st.ntrackr,"ntrackr/I");
+ ntuple_global->Branch("istatr",&ntuple_st.istatr[0],"istatr[500]/I");
+ ntuple_global->Branch("isignr",&ntuple_st.isignr[0],"isignr[500]/I");
+ ntuple_global->Branch("pxr",&ntuple_st.pxr[0],"pxr[500]/F");
+ ntuple_global->Branch("pyr",&ntuple_st.pyr[0],"pyr[500]/F");
+ ntuple_global->Branch("pzr",&ntuple_st.pzr[0],"pzr[500]/F");
+ ntuple_global->Branch("zvr",&ntuple_st.zvr[0],"zvr[500]/F");
+ ntuple_global->Branch("chi2r",&ntuple_st.chi2r[0],"chi2r[500]/F");
+ ntuple_global->Branch("pxv",&ntuple_st.pxv[0],"pxv[500]/F");
+ ntuple_global->Branch("pyv",&ntuple_st.pyv[0],"pyv[500]/F");
+ ntuple_global->Branch("pzv",&ntuple_st.pzv[0],"pzv[500]/F");
+
+ // test aliroot
+
+ new TH1F("h100","particule id du hit geant",20,0.,20.);
+ new TH1F("h101","position en x du hit geant",100,-200.,200.);
+ new TH1F("h102","position en y du hit geant",100,-200.,200.);
+ new TH1F("h103","chambre de tracking concernee",15,0.,14.);
+ new TH1F("h104","moment ptot du hit geant",50,0.,100.);
+ new TH1F("h105","px au vertex",50,0.,20.);
+ new TH1F("h106","py au vertex",50,0.,20.);
+ new TH1F("h107","pz au vertex",50,0.,20.);
+ new TH1F("h108","position zv",50,-15.,15.);
+ new TH1F("h109","position en x du hit reconstruit",100,-300.,300.);
+ new TH1F("h110","position en y du hit reconstruit",100,-300.,300.);
+ new TH1F("h111","delta x ",100,-0.4,0.4);
+ new TH1F("h112","delta y ",100,-0.4,0.4);
+
+ char hname[30];
+ char hname1[30];
+ for (int i=0;i<10;i++) {
+ sprintf(hname,"deltax%d",i);
+ sprintf(hname1,"h12%d",i);
+ new TH1F(hname1,hname ,100,-0.4,0.4);
+ sprintf(hname,"deltay%d",i);
+ sprintf(hname1,"h13%d",i);
+ new TH1F(hname1,hname ,100,-0.4,0.4);
+ }
+ new TH2F("h2000","VAR X st. 5",30,3.0,183.0,100,0.,25.);
+ new TH2F("h2001","VAR Y st. 5",30,3.0,183.0,100,0.,25.);
+
+ new TH2F("h2500","P vs X HHIT",30,3.0,183.0,200,0.,200.);
+ new TH2F("h2501","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.);
+ new TH2F("h2502","P vs X EPH2 st. 5",30,3.0,183.0,100,0.,0.000005);
+ new TH2F("h2503","P vs X EAL2 st. 5",30,3.0,183.0,100,0.,0.01);
+ //new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,0.1);
+ new TH2F("h2505","P vs X EYM2 st. 5",30,3.0,183.0,100,0.,30.);
+
+ new TH2F("h2507","P vs X EPH st. 5",30,3.0,183.0,100,0.,0.003);
+ new TH2F("h2508","P vs X EAL st. 5",30,3.0,183.0,100,0.,0.3);
+ //new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,0.4);
+ new TH2F("h2510","P vs X EYM st. 5",30,3.0,183.0,100,0.,30.);
+
+ new TH2F("h2511","P vs X EPH cut st. 5",30,3.0,183.0,100,0.,0.01);
+ new TH2F("h2512","P vs X EAL cut st. 5",30,3.0,183.0,100,0.,0.3);
+ //new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,0.4);
+ new TH2F("h2514","P vs X EYM cut st. 5",30,3.0,183.0,100,0.,30.);
+ // 4
+ new TH2F("h2400","P vs X HHIT",30,3.0,183.0,200,0.,200.);
+ new TH2F("h2401","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.);
+ new TH2F("h2402","P vs X EPH2 st. 4",30,3.0,183.0,100,0.,0.000005);
+ new TH2F("h2403","P vs X EAL2 st. 4",30,3.0,183.0,100,0.,0.05);
+ //new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,0.1);
+ new TH2F("h2405","P vs X EYM2 st. 4",30,3.0,183.0,100,0.,30.);
+
+ new TH2F("h2407","P vs X EPH st. 4",30,3.0,183.0,100,0.,0.003);
+ new TH2F("h2408","P vs X EAL st. 4",30,3.0,183.0,100,0.,0.3);
+ //new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,0.1);
+ new TH2F("h2410","P vs X EYM st. 4",30,3.0,183.0,100,0.,30.);
+
+ new TH2F("h2411","P vs X EPH cut st. 4",30,3.0,183.0,100,0.,0.01);
+ new TH2F("h2412","P vs X EAL cut st. 4",30,3.0,183.0,100,0.,0.3);
+ //new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,0.1);
+ new TH2F("h2414","P vs X EYM cut st. 4",30,3.0,183.0,100,0.,30.);
+ // 3
+ new TH1F("h2301","P2",30,3.0,183.0);
+ new TH2F("h2302","P2 vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2303","P2 vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.0005);
+ //new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.);
+ new TH2F("h2305","P2 vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.);
+
+ new TH2F("h2307","P vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2308","P vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.);
+ new TH2F("h2310","P vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.);
+
+ new TH2F("h2311","P vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06);
+ new TH2F("h2312","P vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05);
+ //new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.);
+ new TH2F("h2314","P vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.);
+
+ new TH2F("h2315","P2 vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06);
+ new TH2F("h2316","P2 vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05);
+ //new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.);
+ new TH2F("h2318","P2 vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.);
+
+ // 2
+ new TH1F("h2201","P2",30,3.0,183.0);
+ new TH2F("h2202","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2203","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
+ new TH2F("h2205","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.);
+
+ new TH2F("h2207","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2208","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
+ new TH2F("h2210","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.);
+
+ new TH2F("h2211","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05);
+ new TH2F("h2212","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
+ new TH2F("h2214","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.);
+
+ new TH2F("h2215","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05);
+ new TH2F("h2216","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
+ new TH2F("h2218","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.);
+
+ // 1
+ new TH2F("h2102","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2103","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
+ new TH2F("h2105","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.);
+
+ new TH2F("h2107","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2108","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
+ new TH2F("h2110","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.);
+
+ new TH2F("h2111","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1);
+ new TH2F("h2112","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
+ new TH2F("h2114","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.);
+
+ new TH2F("h2115","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1);
+ new TH2F("h2116","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
+ new TH2F("h2118","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.);
+
+ // 2,3,4,5
+ new TH1F("h2701","P2 fit 2",30,3.0,183.0);
+ new TH2F("h2702","P2 vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2703","P2 vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005);
+ // new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.);
+ new TH2F("h2705","P2 vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.);
+
+ new TH2F("h2707","P vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2708","P vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.);
+ new TH2F("h2710","P vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.);
+
+ new TH2F("h2711","P vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07);
+ new TH2F("h2712","P vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.);
+ new TH2F("h2714","P vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.);
+
+ new TH2F("h2715","P2 vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07);
+ new TH2F("h2716","P2 vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.);
+ new TH2F("h2718","P2 vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.);
+
+ // 1,3,4,5
+ new TH1F("h2801","P2 fit 1",30,3.0,183.0);
+ new TH2F("h2802","P2 vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2803","P2 vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.);
+ new TH2F("h2805","P2 vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.);
+
+ new TH2F("h2807","P vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h2808","P vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005);
+ //new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.);
+ new TH2F("h2810","P vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.);
+
+ new TH2F("h2811","P vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05);
+ new TH2F("h2812","P vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.);
+ new TH2F("h2814","P vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.);
+
+ new TH2F("h2815","P2 vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05);
+ new TH2F("h2816","P2 vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2);
+ //new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
+ new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.);
+ new TH2F("h2818","P2 vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.);
+ // fin de test
+
+ new TH1F("h500","Acceptance en H st. 4",500,0.,500.);
+ new TH1F("h600","Acceptance en H st. 5",500,0.,500.);
+ new TH1F("h700","X vertex track found",200,-10.,10.);
+ new TH1F("h701","Y vertex track found",200,-10.,10.);
+ new TH1F("h800","Rap. muon gen.",100,0.,5.);
+ new TH1F("h801","Rap. muon gen. recons.",100,0.,5.);
+ new TH1F("h802","Rap. muon gen. ghost ",100,0.,5.);
+ new TH1F("h900","Pt muon gen.",100,0.,20.);
+ new TH1F("h901","Pt muon gen. recons.",100,0.,20.);
+ new TH1F("h902","Pt muon gen. ghost",100,0.,20.);
+ new TH1F("h910","phi muon gen.",100,-10.,10.);
+ new TH1F("h911","phi muon gen. recons.",100,-10.,10.);
+ new TH1F("h912","phi muon gen. ghost",100,-10.,10.);
+ new TH2F("h1001","Y VS X hit st. 1",300,-300.,300.,300,-300.,300.);
+ new TH2F("h1002","Y VS X hit st. 2",300,-300.,300.,300,-300.,300.);
+ new TH2F("h1003","Y VS X hit st. 3",300,-300.,300.,300,-300.,300.);
+ new TH2F("h1004","Y VS X hit st. 4",300,-300.,300.,300,-300.,300.);
+ new TH2F("h1005","Y VS X hit st. 5",300,-300.,300.,300,-300.,300.);
+ // Histos variance dans 4
+ new TH2F("h11","VAR X st. 4",30,3.0,183.0,100,0.,2.);
+ new TH2F("h12","VAR Y st. 4",30,3.0,183.0,100,0.,600.);
+ new TH2F("h13","VAR PHI st. 4",30,3.0,183.0,100,0.,0.0001);
+ new TH2F("h14","VAR ALM st. 4",30,3.0,183.0,100,0.,0.05);
+ new TH1F("h15","P",30,3.0,183.0);
+ new TH1F("h411","VAR X st. 4",100,-1.42,1.42);
+ new TH1F("h412","VAR Y st. 4",100,-25.,25.);
+ new TH1F("h413","VAR PHI st. 4",100,-0.01,0.01);
+ new TH1F("h414","VAR ALM st. 4",100,-0.23,0.23);
+ // histo2
+ new TH2F("h211","histo2-VAR X st. 4",30,3.0,183.0,100,0.,2.);
+ new TH2F("h212","histo2-VAR Y st. 4",30,3.0,183.0,100,0.,600.);
+ new TH1F("h213","histo2-VAR X st. 4",100,-1.42,1.42);
+ new TH1F("h214","histo2-VAR Y st. 4",100,-25.,25.);
+ new TH1F("h215","histo2-P",30,3.0,183.0);
+
+ // Histos variance dans 2
+ new TH2F("h21","VAR X st. 2",30,3.0,183.0,100,0.,3.);
+ new TH2F("h22","VAR Y st. 2",30,3.0,183.0,100,0.,7.);
+ new TH2F("h23","VAR PHI st. 2",30,3.0,183.0,100,0.,0.006);
+ new TH2F("h24","VAR ALM st. 2",30,3.0,183.0,100,0.,0.005);
+ new TH1F("h25","P",30,3.0,183.0);
+ new TH1F("h421","VAR X st. 2",100,-1.72,1.72);
+ new TH1F("h422","VAR Y st. 2",100,-2.7,2.7);
+ new TH1F("h423","VAR PHI st. 2",100,-0.08,0.08);
+ new TH1F("h424","VAR ALM st. 2",100,-0.072,0.072);
+ // histo2
+ new TH2F("h221","histo2-VAR X st. 2",30,3.0,183.0,100,0.,3.);
+ new TH2F("h222","histo2-VAR Y st. 2",30,3.0,183.0,100,0.,7.);
+ new TH1F("h223","histo2-VAR X st. 2",100,-1.72,1.72);
+ new TH1F("h224","histo2-VAR Y st. 2",100,-2.7,2.7);
+ new TH1F("h225","histo2-P",30,3.0,183.0);
+
+ // Histos variance dans 1
+ new TH2F("h31","VAR X st. 1",30,3.0,183.0,100,0.,2.);
+ new TH2F("h32","VAR Y st. 1",30,3.0,183.0,100,0.,0.5);
+ new TH2F("h33","VAR PHI st. 1",30,3.0,183.0,100,0.,0.006);
+ new TH2F("h34","VAR ALM st. 1",30,3.0,183.0,100,0.,0.005);
+ new TH1F("h35","P",30,3.0,183.0);
+ new TH1F("h431","VAR X st. 1",100,-1.42,1.42);
+ new TH1F("h432","VAR Y st. 1",100,-0.72,0.72);
+ new TH1F("h433","VAR PHI st. 1",100,-0.08,0.08);
+ new TH1F("h434","VAR ALM st. 1",100,-0.072,0.072);
+ // Histos variance dans 1
+ new TH2F("h41","VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.);
+ new TH2F("h42","VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.);
+ new TH2F("h43","VAR PHI st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005);
+ new TH2F("h44","VAR ALM st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005);
+ new TH1F("h45","P",30,3.0,183.0);
+ new TH1F("h441","VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.);
+ new TH1F("h442","VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5);
+ new TH1F("h443","VAR PHI st. 1 fit 5,4,3,2,V",100,-0.072,0.072);
+ new TH1F("h444","VAR ALM st. 1 fit 5,4,3,2,V",100,-0.072,0.072);
+ // histo2
+ new TH2F("h241","histo2-VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.);
+ new TH2F("h242","histo2-VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.);
+ new TH1F("h243","histo2-VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.);
+ new TH1F("h244","histo2-VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5);
+ new TH1F("h245","histo2-P",30,3.0,183.0);
+
+ // Histos variance dans 2
+ new TH2F("h51","VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5);
+ new TH2F("h52","VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.);
+ new TH2F("h53","VAR PHI st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.005);
+ new TH2F("h54","VAR ALM st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.01);
+ new TH1F("h55","P",30,3.0,183.0);
+ new TH1F("h451","VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72);
+ new TH1F("h452","VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42);
+ new TH1F("h453","VAR PHI st. 2 fit 5,4,3,1,V",100,-0.072,0.072);
+ new TH1F("h454","VAR ALM st. 2 fit 5,4,3,1,V",100,-0.1,0.1);
+ new TH1F("h999","PTOT",30,3.0,183.0);
+ // histo2
+ new TH2F("h251","histo2-VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5);
+ new TH2F("h252","histo2-VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.);
+ new TH1F("h253","histo2-VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72);
+ new TH1F("h254","histo2-VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42);
+ new TH1F("h255","histo2-P",30,3.0,183.0);
+ // Histos variance dans 3
+ new TH2F("h61","VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.);
+ new TH2F("h62","VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.);
+ new TH2F("h63","VAR PHI st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006);
+ new TH2F("h64","VAR ALM st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006);
+ new TH1F("h65","P",30,3.0,183.0);
+ new TH1F("h461","VAR X st. 3 fit 4,5,V",100,-2.25,2.25);
+ new TH1F("h462","VAR Y st. 3 fit 4,5,V",100,-1.42,1.42);
+ new TH1F("h463","VAR PHI st. 3 fit 4,5,V",100,-0.024,0.024);
+ new TH1F("h464","VAR ALM st. 3 fit 4,5,V",100,-0.024,0.024);
+ // histo2
+ new TH2F("h261","histo2-VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.);
+ new TH2F("h262","histo2-VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.);
+ new TH1F("h263","histo2-VAR X st. 3 fit 4,5,V",100,-2.25,2.25);
+ new TH1F("h264","histo2-VAR Y st. 3 fit 4,5,V",100,-1.42,1.42);
+ new TH1F("h265","Phisto2-",30,3.0,183.0);
+ // Histos dx,dy distribution between chambers inside stations
+ new TH1F("h71","DX in st. ID-70",100,-5.,5.);
+ new TH1F("h81","DY in st. ID-80",100,-5.,5.);
+ new TH1F("h72","DX in st. ID-70",100,-5.,5.);
+ new TH1F("h82","DY in st. ID-80",100,-5.,5.);
+ new TH1F("h73","DX in st. ID-70",100,-5.,5.);
+ new TH1F("h83","DY in st. ID-80",100,-5.,5.);
+ new TH1F("h74","DX in st. ID-70",100,-5.,5.);
+ new TH1F("h84","DY in st. ID-80",100,-5.,5.);
+ new TH1F("h75","DX in st. ID-70",100,-5.,5.);
+ new TH1F("h85","DY in st. ID-80",100,-5.,5.);
+}
+
+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)
+{
+ //
+ // fill the ntuple
+ ntuple_st.ievr = ievr;
+ ntuple_st.ntrackr = ntrackr;
+ for (Int_t i=0; i<500; i++) {
+ ntuple_st.istatr[i] = istatr[i];
+ ntuple_st.isignr[i] = isignr[i];
+ ntuple_st.pxr[i] = pxr[i];
+ ntuple_st.pyr[i] = pyr[i];
+ ntuple_st.pzr[i] = pzr[i];
+ ntuple_st.zvr[i] = zvr[i];
+ ntuple_st.chi2r[i] = chi2r[i];
+ ntuple_st.pxv[i] = pxv[i];
+ ntuple_st.pyv[i] = pyv[i];
+ ntuple_st.pzv[i] = pzv[i];
+ }
+ ntuple_global->Fill();
+}
+
+void hist_closed()
+{
+ //
+ // write histos and ntuple to "reconst.root" file
+ hfile_global->Write();
+}
+
+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)
+{
+ //
+ // introduce aliroot variables in fortran common
+ // tracking study from geant hits
+ //
+
+ AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
+
+ // TTree *TK = gAlice->TreeK();
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks = (Int_t)TH->GetEntries();
+ cout<<"ntrack="<<ntracks<<endl;
+
+ Int_t maxidg = 0;
+ Int_t nres=0;
+
+//
+// Loop over tracks
+//
+
+ for (Int_t track=0; track<ntracks;track++) {
+ gAlice->ResetHits();
+ TH->GetEvent(track);
+
+ if (MUON) {
+//
+// Loop over hits
+//
+ for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
+ mHit;
+ mHit=(AliMUONhit*)MUON->NextHit())
+ {
+ if (maxidg<=20000) {
+
+ if (mHit->fChamber > 10) continue;
+ TClonesArray *fPartArray = gAlice->Particles();
+ TParticle *Part;
+ Int_t ftrack = mHit->fTrack;
+ Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
+
+ if (id==kMuonPlus||id==kMuonMinus) {
+
+ // inversion de x et y car le champ est inverse dans le programme tracking
+ xtrg[maxidg] = 0;
+ ytrg[maxidg] = 0;
+ xgeant[maxidg] = mHit->fY; // x-pos of hit
+ ygeant[maxidg] = mHit->fX; // y-pos of hit
+ clsize1[maxidg] = 0; // cluster size on 1-st cathode
+ clsize2[maxidg] = 0; // cluster size on 2-nd cathode
+ cx[maxidg] = mHit->fCyHit; // Px/P of hit
+ cy[maxidg] = mHit->fCxHit; // Py/P of hit
+ cz[maxidg] = mHit->fCzHit; // Pz/P of hit
+ izch[maxidg] = mHit->fChamber;
+ /*
+ Int_t pdgtype = Int_t(mHit->fParticle); // particle number
+ itypg[maxidg] = gMC->IdFromPDG(pdgtype);
+
+ */
+ if (id==kMuonPlus) itypg[maxidg] = 5;
+ else itypg[maxidg] = 6;
+
+ ptotg[maxidg] = mHit->fPTot; // P of hit
+
+ Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
+ Float_t thet = Part->Theta();
+ thet = thet*180./3.1416;
+
+ Int_t iparent = Part->GetFirstMother();
+ if (iparent >= 0) {
+ Int_t ip;
+ while(1) {
+ ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother();
+ if (ip < 0) {
+ break;
+ } else {
+ iparent = ip;
+ }
+ }
+ }
+ //printf("iparent - %d\n",iparent);
+ Int_t id1 = ftrack; // numero de la particule generee au vertex
+ Int_t idum = track+1;
+ Int_t id2 = ((TParticle*) fPartArray->UncheckedAt(iparent))->GetPdgCode();
+
+ if (id2==443) id2=114;
+ else id2=116;
+
+ if (id2==116) {
+ nres++;
+ }
+ //printf("id2 %d\n",id2);
+ idg[maxidg] = 30000*id1+10000*idum+id2;
+
+ pvert1g[maxidg] = Part->Py(); // Px vertex
+ pvert2g[maxidg] = Part->Px(); // Py vertex
+ pvert3g[maxidg] = Part->Pz(); // Pz vertex
+ zvertg[maxidg] = Part->Vz(); // z vertex
+ maxidg ++;
+
+ }
+ }
+ } // hit loop
+ } // if MUON
+ } // track loop first file
+
+ if (TrH1 && fHits2 ) { // if background file
+ ntracks =(Int_t)TrH1->GetEntries();
+ printf("Trackf_read - 2-nd file - ntracks %d\n",ntracks);
+
+ // Loop over tracks
+ for (Int_t track=0; track<ntracks; track++) {
+
+ if (fHits2) fHits2->Clear();
+ TrH1->GetEvent(track);
+
+ // Loop over hits
+ for (int i=0;i<fHits2->GetEntriesFast();i++)
+ {
+ AliMUONhit *mHit=(AliMUONhit*) (*fHits2)[i];
+
+ if (mHit->fChamber > 10) continue;
+
+ if (maxidg<=20000) {
+
+ // inversion de x et y car le champ est inverse dans le programme tracking !!!!
+ xtrg[maxidg] = 0; // only for reconstructed point
+ ytrg[maxidg] = 0; // only for reconstructed point
+ xgeant[maxidg] = mHit->fY; // x-pos of hit
+ ygeant[maxidg] = mHit->fX; // y-pos of hit
+ clsize1[maxidg] = 0; // cluster size on 1-st cathode
+ clsize2[maxidg] = 0; // cluster size on 2-nd cathode
+ cx[maxidg] = mHit->fCyHit; // Px/P of hit
+ cy[maxidg] = mHit->fCxHit; // Py/P of hit
+ cz[maxidg] = mHit->fCzHit; // Pz/P of hit
+ izch[maxidg] = mHit->fChamber; // chamber number
+ ptotg[maxidg] = mHit->fPTot; // P of hit
+
+ Int_t ftrack = mHit->fTrack;
+ Int_t id1 = ftrack; // track number
+ Int_t idum = track+1;
+
+ TClonesArray *fPartArray = fParticles2;
+ TParticle *Part;
+ Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
+ Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
+ if (id==kMuonPlus||id==kMuonMinus) {
+ if (id==kMuonPlus) itypg[maxidg] = 5;
+ else itypg[maxidg] = 6;
+ } else itypg[maxidg]=0;
+
+ Int_t id2=0; // set parent to 0 for background !!
+ idg[maxidg] = 30000*id1+10000*idum+id2;
+
+ pvert1g[maxidg] = Part->Py(); // Px vertex
+ pvert2g[maxidg] = Part->Px(); // Py vertex
+ pvert3g[maxidg] = Part->Pz(); // Pz vertex
+ zvertg[maxidg] = Part->Vz(); // z vertex
+
+ maxidg ++;
+
+ } // check limits (maxidg)
+ } // hit loop
+ } // track loop
+ } // if TrH1
+
+ ievr = nev;
+ nhittot1 = maxidg ;
+ cout<<"nhittot1="<<nhittot1<<endl;
+
+ static Int_t nbres=0;
+ if (nres>=19) nbres++;
+ printf("nres, nbres %d %d \n",nres,nbres);
+
+ hfile_global->cd();
+
+}
+
+
+
+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)
+
+{
+ //
+ // introduce aliroot variables in fortran common
+ // tracking study from reconstructed points
+ //
+ AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
+
+ cout<<"numero de l'evenement "<<nev<<endl;
+
+ MUON->GetTreeC(nev);
+ TTree *TC=MUON->TreeC();
+ TC->GetEntries();
+
+ Int_t maxidg = 0;
+ Int_t nres=0;
+ Int_t nncor=0;
+ static Int_t nuncor=0;
+ static Int_t nbadcor=0;
+ AliMUONRawCluster * mRaw;
+ AliMUONRawCluster * mRaw1;
+ TTree *TH = gAlice->TreeH();
+
+ Int_t ihit;
+ Int_t mult1, mult2;
+ if (MUON) {
+ for (Int_t ich=0;ich<10;ich++) {
+ TClonesArray *MUONcorrel = MUON->CathCorrelAddress(ich);
+ MUON->ResetCorrelation();
+ TC->GetEvent();
+ Int_t ncor = (Int_t)MUONcorrel->GetEntries();
+ if (ncor>=2) nncor++;
+ if (!ncor) continue;
+
+ // Loop over correlated clusters
+ for (Int_t icor=0;icor<ncor;icor++) {
+ AliMUONcorrelation * mCor = (AliMUONcorrelation*)MUONcorrel->UncheckedAt(icor);
+
+ Int_t flag=0; // = 1 if no information in the second cathode
+ Int_t index = mCor->fCorrelIndex[0]; // for the second cathode
+ if (index >= 0) {
+ Int_t index1 = mCor->fCorrelIndex[3]; // for the 1-st cathode
+ mRaw1 = MUON->RawCluster(ich,1,index1);
+ mult1=mRaw1->fMultiplicity;
+ mRaw = MUON->RawCluster(ich,2,index);
+ mult2=mRaw->fMultiplicity;
+ } else {
+ index = mCor->fCorrelIndex[3];
+ mRaw = MUON->RawCluster(ich,1,index);
+ mult1=mRaw->fMultiplicity;
+ mult2=0;
+ flag=1;
+ nuncor++;
+ }
+ if (!mRaw) continue;
+
+ Int_t ftrack1 = mRaw->fTracks[1]; // qui doit etre le meme pour
+ // la cathode 1 et 2
+ ihit= mRaw->fTracks[0];
+ //printf("icor, ftrack1 ihit %d %d %d\n",icor,ftrack1,ihit);
+
+ if (mRaw->fClusterType == 0 ) {
+
+ if (maxidg<=20000) {
+ if (flag == 0) {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[0];
+ Int_t index1 = mCor->fCorrelIndex[3];
+ mRaw1 = MUON->RawCluster(ich,1,index1);
+ if (mRaw1->fClusterType==1 || mRaw1->fClusterType==2) {
+ Float_t xclust=mCor->fX[3];
+ Float_t yclust=mCor->fY[3];
+ AliMUONchamber *iChamber=&(MUON->Chamber(ich));
+ AliMUONsegmentation *seg = iChamber->GetSegmentationModel(1);
+ Int_t ix,iy;
+ seg->GetPadIxy(xclust,yclust,ix,iy);
+ Int_t isec=seg->Sector(ix,iy);
+ printf("nev, CORRELATION with pure background in chamber sector %d %d %d !!!!!!!!!!!!!!!!!!!!!\n",nev,ich+1,isec);
+ nbadcor++;
+
+ } // end if cluster type on cathode 1
+ }else {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[3];
+ } // if iflag
+ izch[maxidg] = ich+1;
+ xgeant[maxidg] = 0;
+ ygeant[maxidg] = 0;
+ clsize1[maxidg] = mult1;
+ clsize2[maxidg] = mult2;
+
+ cx[maxidg] = 0; // Px/P of hit
+ cy[maxidg] = 0; // Py/P of hit
+ cz[maxidg] = 0; // Pz/P of hit
+ itypg[maxidg] = 0; // particle number
+ ptotg[maxidg] = 0; // P of hit
+ idg[maxidg] = 0;
+ pvert1g[maxidg] = 0; // Px vertex
+ pvert2g[maxidg] = 0; // Py vertex
+ pvert3g[maxidg] = 0; // Pz vertex
+ zvertg[maxidg] = 0; // z vertex
+ maxidg++;
+
+ }// fin maxidg
+
+ } else if (mRaw->fClusterType ==1 && ftrack1 < 0) // background + resonance
+ {
+ nres++;
+ // get indexmap and loop over digits to find the signal
+ Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
+ gAlice->ResetDigits();
+ if (flag==0) {
+ //gAlice->TreeD()->GetEvent(2); // cathode 2
+ gAlice->TreeD()->GetEvent(nent-1); // cathode 2
+ } else {
+ //gAlice->TreeD()->GetEvent(1); // cathode 1
+ gAlice->TreeD()->GetEvent(nent-2); // cathode 1
+ }
+
+ TClonesArray *MUONdigits = MUON->DigitsAddress(ich);
+ Int_t mul=mRaw->fMultiplicity;
+ Int_t trsign;
+ for (int i=0;i<mul;i++) {
+ Int_t idx=mRaw->fIndexMap[i];
+ AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx);
+ trsign=dig->fTracks[0];
+ ihit=dig->fHit-1;
+ if (trsign > 0 && ihit >= 0) break;
+
+ } // loop over indexmap
+
+ //printf("trsign, ihit %d %d\n",trsign,ihit);
+ //printf("signal+background : trsign %d\n",trsign);
+
+ if (trsign < 0 || ihit < 0) { // no signal muon was found
+
+ if (maxidg<=20000) {
+ if (flag == 0) {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[0];
+ }else {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[3];
+ }
+
+ izch[maxidg] = ich+1;
+
+ // initialisation of informations which
+ // can't be reached for background
+
+ xgeant[maxidg] = 0; // only for resonances
+ ygeant[maxidg] = 0; // only for resonances
+ clsize1[maxidg] = mult1;
+ clsize2[maxidg] = mult2;
+
+ cx[maxidg] = 0; // Px/P of hit
+ cy[maxidg] = 0; // Py/P of hit
+ cz[maxidg] = 0; // Pz/P of hit
+ itypg[maxidg] = 0; // particle number
+ ptotg[maxidg] = 0; // P of hit
+ idg[maxidg] = 0;
+ pvert1g[maxidg] = 0; // Px vertex
+ pvert2g[maxidg] = 0; // Py vertex
+ pvert3g[maxidg] = 0; // Pz vertex
+ zvertg[maxidg] = 0;
+ maxidg++;
+
+ }// fin maxidg
+ } else { // signal muon - retrieve info
+ //printf("inside trsign, ihit %d %d\n",trsign,ihit);
+ if (maxidg<=20000) {
+ if (flag == 0) {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[0];
+ }else {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[3];
+ }
+ izch[maxidg] = ich+1;
+ clsize1[maxidg] = mult1;
+ clsize2[maxidg] = mult2;
+
+ // initialise and set to the correct values
+ // if signal muons
+
+ xgeant[maxidg] = 0; // only for resonances
+ ygeant[maxidg] = 0; // only for resonances
+
+ cx[maxidg] = 0; // Px/P of hit
+ cy[maxidg] = 0; // Py/P of hit
+ cz[maxidg] = 0; // Pz/P of hit
+ itypg[maxidg] = 0; // particle number
+ ptotg[maxidg] = 0; // P of hit
+ idg[maxidg] = 0;
+ pvert1g[maxidg] = 0; // Px vertex
+ pvert2g[maxidg] = 0; // Py vertex
+ pvert3g[maxidg] = 0; // Pz vertex
+ zvertg[maxidg] = 0;
+ // try to retrieve info about signal muons
+ gAlice->ResetHits();
+ TH->GetEvent(trsign);
+
+ TClonesArray *MUONhits = MUON->Hits();
+ AliMUONhit *mHit= (AliMUONhit*)MUONhits->
+ UncheckedAt(ihit);
+ TClonesArray *fPartArray = gAlice->Particles();
+ TParticle *Part;
+ Int_t nch=mHit->fChamber-1;
+ //printf("sig+bgr ich, nch %d %d \n",ich,nch);
+ if (nch==ich) {
+ Int_t ftrack = mHit->fTrack;
+ Int_t id = ((TParticle*) fPartArray->
+ UncheckedAt(ftrack))->GetPdgCode();
+ if (id==kMuonPlus||id==kMuonMinus) {
+ xgeant[maxidg] = (Double_t) mHit->fY;
+ ygeant[maxidg] = (Double_t) mHit->fX;
+ cx[maxidg] = (Double_t) mHit->fCyHit;
+ cy[maxidg] = (Double_t) mHit->fCxHit;
+ cz[maxidg] = (Double_t) mHit->fCzHit;
+
+ if (id==kMuonPlus) {
+ itypg[maxidg] = 5;
+ } else if (id==kMuonMinus) {
+ itypg[maxidg] = 6;
+ } else itypg[maxidg] = 0;
+
+ ptotg[maxidg] = (Double_t) mHit->fPTot;
+ Part = (TParticle*) fPartArray->
+ UncheckedAt(ftrack);
+ Int_t iparent = Part->GetFirstMother();
+ Int_t id2;
+ id2 = ((TParticle*) fPartArray->
+ UncheckedAt(ftrack))->GetPdgCode();
+
+ if (iparent >= 0) {
+ Int_t ip;
+ while(1) {
+ ip=((TParticle*) fPartArray->
+ UncheckedAt(iparent))->GetFirstMother();
+ if (ip < 0) {
+ id2 = ((TParticle*) fPartArray->
+ UncheckedAt(iparent))->GetPdgCode();
+ break;
+ } else {
+ iparent = ip;
+ id2 = ((TParticle*) fPartArray->
+ UncheckedAt(iparent))->GetPdgCode();
+ } // ip<0
+ } // while
+ }// iparent
+ Int_t id1 = ftrack;
+ Int_t idum = trsign+1;
+
+ if (id2==443 || id2==553) {
+ nres++;
+ if (id2==443) id2=114;
+ else id2=116;
+ }
+
+ idg[maxidg] = 30000*id1+10000*idum+id2;
+ pvert1g[maxidg] = (Double_t) Part->Py();
+ pvert2g[maxidg] = (Double_t) Part->Px();
+ pvert3g[maxidg] = (Double_t) Part->Pz();
+ zvertg[maxidg] = (Double_t) Part->Vz();
+ } //if muon
+ } //if nch
+ maxidg++;
+ } // check limits
+ } // sign+bgr, highest bgr
+ }
+ //pure resonance or mixed cluster with the highest
+ //contribution coming from resonance
+ if (mRaw->fClusterType >= 1 && ftrack1>=0)
+ {
+ if (maxidg<=20000) {
+ if (flag == 0) {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[0];
+ }else {
+ xtrg[maxidg] = (Double_t) mCor->fY[3];
+ ytrg[maxidg] = (Double_t) mCor->fX[3];
+ }
+ clsize1[maxidg] = mult1;
+ clsize2[maxidg] = mult2;
+ izch[maxidg] = ich+1;
+
+ Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
+ gAlice->ResetDigits();
+ if (flag==0) {
+ //gAlice->TreeD()->GetEvent(2); // cathode 2
+ gAlice->TreeD()->GetEvent(nent-1); // cathode 2
+ } else {
+ //gAlice->TreeD()->GetEvent(1); // cathode 1
+ gAlice->TreeD()->GetEvent(nent-2); // cathode 1
+ }
+
+ TClonesArray *MUONdigits = MUON->DigitsAddress(ich);
+ Int_t mul=mRaw->fMultiplicity;
+ for (int i=0;i<mul;i++) {
+ Int_t idx=mRaw->fIndexMap[i];
+ AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx);
+ ihit=dig->fHit-1;
+ if (ihit >= 0) break;
+
+ } // loop over indexmap
+ //printf("fClusterType, ihit %d %d \n",mRaw->fClusterType,ihit);
+ if (ihit < 0) {
+ xgeant[maxidg] = 0; // only for resonances
+ ygeant[maxidg] = 0; // only for resonances
+
+ cx[maxidg] = 0; // Px/P of hit
+ cy[maxidg] = 0; // Py/P of hit
+ cz[maxidg] = 0; // Pz/P of hit
+ itypg[maxidg] = 0; // particle number
+ ptotg[maxidg] = 0; // P of hit
+ idg[maxidg] = 0;
+ pvert1g[maxidg] = 0; // Px vertex
+ pvert2g[maxidg] = 0; // Py vertex
+ pvert3g[maxidg] = 0; // Pz vertex
+ zvertg[maxidg] = 0;
+ } else {
+ gAlice->ResetHits();
+ TH->GetEvent(ftrack1);
+ TClonesArray *MUONhits = MUON->Hits();
+ AliMUONhit *mHit= (AliMUONhit*)MUONhits->
+ UncheckedAt(ihit);
+ TClonesArray *fPartArray = gAlice->Particles();
+ TParticle *Part;
+ Int_t nch=mHit->fChamber-1;
+ //printf("signal ich, nch %d %d \n",ich,nch);
+ if (nch==ich) {
+ Int_t ftrack = mHit->fTrack;
+ Int_t id = ((TParticle*) fPartArray->
+ UncheckedAt(ftrack))->GetPdgCode();
+ //printf("id %d \n",id);
+ if (id==kMuonPlus||id==kMuonMinus) {
+ xgeant[maxidg] = (Double_t) mHit->fY;
+ ygeant[maxidg] = (Double_t) mHit->fX;
+ cx[maxidg] = (Double_t) mHit->fCyHit;
+ cy[maxidg] = (Double_t) mHit->fCxHit;
+ cz[maxidg] = (Double_t) mHit->fCzHit;
+
+ if (id==kMuonPlus) {
+ itypg[maxidg] = 5;
+ } else if (id==kMuonMinus) {
+ itypg[maxidg] = 6;
+ } else itypg[maxidg] = 0;
+
+ ptotg[maxidg] = (Double_t) mHit->fPTot;
+ Part = (TParticle*) fPartArray->
+ UncheckedAt(ftrack);
+ Int_t iparent = Part->GetFirstMother();
+ Int_t id2;
+ id2 = ((TParticle*) fPartArray->
+ UncheckedAt(ftrack))->GetPdgCode();
+
+ if (iparent >= 0) {
+ Int_t ip;
+ while(1) {
+ ip=((TParticle*) fPartArray->
+ UncheckedAt(iparent))->GetFirstMother();
+ if (ip < 0) {
+ id2 = ((TParticle*) fPartArray->
+ UncheckedAt(iparent))->GetPdgCode();
+ break;
+ } else {
+ iparent = ip;
+ id2 = ((TParticle*) fPartArray->
+ UncheckedAt(iparent))->GetPdgCode();
+ } // ip<0
+ } // while
+ }// iparent
+ Int_t id1 = ftrack;
+ Int_t idum = ftrack1+1;
+
+ if (id2==443 || id2==553) {
+ nres++;
+ if (id2==443) id2=114;
+ else id2=116;
+ }
+ // printf("id2 %d\n",id2);
+ idg[maxidg] = 30000*id1+10000*idum+id2;
+ pvert1g[maxidg] = (Double_t) Part->Py();
+ pvert2g[maxidg] = (Double_t) Part->Px();
+ pvert3g[maxidg] = (Double_t) Part->Pz();
+ zvertg[maxidg] = (Double_t) Part->Vz();
+ } //if muon
+ } //if nch
+ } // ihit
+ maxidg++;
+ } // check limits
+ } // if cluster type
+ } // icor loop
+ } // ich loop
+ }// if MUON
+
+
+ ievr = nev;
+ cout<<"evenement "<<ievr<<endl;
+ nhittot1 = maxidg ;
+ cout<<"nhittot1="<<nhittot1<<endl;
+
+ static Int_t nbres=0;
+ static Int_t nbcor=0;
+ if (nres>=19) nbres++;
+ printf("nres ,nncor - %d %d\n",nres,nncor);
+ printf("nbres - %d\n",nbres);
+ if (nncor>=20) nbcor++;
+ printf("nbcor - %d\n",nbcor);
+ printf("nuncor - %d\n",nuncor);
+ printf("nbadcor - %d\n",nbadcor);
+
+ TC->Reset();
+
+ hfile_global->cd();
+
+}
+
+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)
+{
+ //
+ // Fit a track candidate with the following input parameters:
+ // INPUT : IVERTEX : vertex flag, if IVERTEX=1 (XVERT,YVERT) are free paramaters
+ // if IVERTEX=1 (XVERT,YVERT)=(0.,0.)
+ // PEST(5) : starting value of parameters (minuit)
+ // PSTEP(5) : step size for parameters (minuit)
+ // OUTPUT : PXZINV,TPHI,TALAM,XVERT,YVERT : fitted value of the parameters
+
+ static Double_t arglist[10];
+ static Double_t c[5] = {0.4, 0.45, 0.45, 90., 90.};
+ static Double_t b1, b2, epxz, efi, exs, exvert, eyvert;
+ TString chname;
+ Int_t ierflg = 0;
+
+ TMinuit *gMinuit = new TMinuit(5);
+ gMinuit->mninit(5,10,7);
+ gMinuit->SetFCN(fcnf); // constant m.f.
+
+ arglist[0] = -1;
+
+ gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
+ // gMinuit->mnseti('track fitting');
+
+ gMinuit->mnparm(0, "invmom", pest[0], pstep[0], -c[0], c[0], ierflg);
+ gMinuit->mnparm(1, "azimuth", pest[1], pstep[1], -c[1], c[1], ierflg);
+ gMinuit->mnparm(2, "deep", pest[2], pstep[2], -c[2], c[2], ierflg);
+ if (ivertex==1) {
+ gMinuit->mnparm(3, "x ", pest[3], pstep[3], -c[3], c[3], ierflg);
+ gMinuit->mnparm(4, "y ", pest[4], pstep[4], -c[4], c[4], ierflg);
+ }
+
+ gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg);
+ gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg);
+ gMinuit->mnexcm("EXIT" , arglist, 0, ierflg);
+
+ gMinuit->mnpout(0, chname, pxzinv, epxz, b1, b2, ierflg);
+ gMinuit->mnpout(1, chname, tphi, efi, b1, b2, ierflg);
+ gMinuit->mnpout(2, chname, talam, exs, b1, b2, ierflg);
+ if (ivertex==1) {
+ gMinuit->mnpout(3, chname, xvert, exvert, b1, b2, ierflg);
+ gMinuit->mnpout(4, chname, yvert, eyvert, b1, b2, ierflg);
+ }
+
+ delete gMinuit;
+
+}
+
+void fcnf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *pest, Int_t iflag)
+{
+ //
+ // function called by trackf_fit
+ Int_t futil = 0;
+ fcn(npar,grad,fval,pest,iflag,futil);
+}
+
+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)
+{
+ //
+ // minuit fits for tracking finding
+
+ static Double_t arglist[10];
+ static Double_t c1[5] = {0.001, 0.001, 0.001, 1., 1.};
+ static Double_t c2[5] = {0.5, 0.5, 0.5, 120., 120.};
+ static Double_t emat[9];
+ static Double_t b1, b2;
+ Double_t fmin, fedm, errdef;
+ Int_t npari, nparx, istat;
+
+ TString chname;
+ Int_t ierflg = 0;
+
+ TMinuit *gMinuit = new TMinuit(5);
+ gMinuit->mninit(5,10,7);
+ gMinuit->SetFCN(fcnfitf);
+
+ arglist[0] = -1.;
+ gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
+
+ // gMinuit->mnseti('track fitting');
+
+ gMinuit->mnparm(0,"invmom", pxzinv, c1[0], -c2[0], c2[0], ierflg); // 0.003, 0.5
+ gMinuit->mnparm(1,"azimuth ", fis, c1[1], -c2[1], c2[1], ierflg);
+ gMinuit->mnparm(2,"deep ", alams, c1[2], -c2[2], c2[2], ierflg);
+ gMinuit->mnparm(3,"xvert", xvert, c1[3], -c2[3], c2[3], ierflg);
+ gMinuit->mnparm(4,"yvert", yvert, c1[4], -c2[4], c2[4], ierflg);
+
+ gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg);
+ arglist[0] = 2.;
+ gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg);
+ gMinuit->mnexcm("EXIT", arglist, 0, ierflg);
+
+ gMinuit->mnpout(0, chname, pxzinvf, epxzinv, b1, b2, ierflg);
+ gMinuit->mnpout(1, chname, fif, efi, b1, b2, ierflg);
+ gMinuit->mnpout(2, chname, alf, exs, b1, b2, ierflg);
+ gMinuit->mnpout(3, chname, xvertf, exvert, b1, b2, ierflg);
+ gMinuit->mnpout(4, chname, yvertf, eyvert, b1, b2, ierflg);
+
+ gMinuit->mnemat(emat, 3);
+ gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
+
+ delete gMinuit;
+}
+
+void fcnfitf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *xval, Int_t iflag)
+{
+ //
+ // function called by prec_fit
+ Int_t futil = 0;
+ fcnfit(npar,grad,fval,xval,iflag,futil);
+}
+
+///////////////////// fin modifs perso //////////////////////
+
ClassImp(AliMUONcluster)
//___________________________________________
//
// Creates a MUON digit object to be updated
//
- fPadX = digits[0];
- fPadY = digits[1];
- fSignal = digits[2];
+ fPadX = digits[0];
+ fPadY = digits[1];
+ fSignal = digits[2];
+ fPhysics = digits[3];
+ fHit = digits[4];
}
//_____________________________________________________________________________
fPadX = digits[0];
fPadY = digits[1];
fSignal = digits[2];
+ fPhysics = digits[3];
+ fHit = digits[4];
for(Int_t i=0; i<10; i++) {
fTcharges[i] = charges[i];
fTracks[i] = tracks[i];
}
}
+AliMUONdigit::~AliMUONdigit()
+{
+
+}
+
ClassImp(AliMUONlist)
//____________________________________________________________________________
- AliMUONlist::AliMUONlist(Int_t rpad, Int_t *digits):
+ AliMUONlist::AliMUONlist(Int_t ich, Int_t *digits):
AliMUONdigit(digits)
{
//
// Creates a MUON digit list object
//
- fRpad = rpad;
+ fChamber = ich;
fTrackList = new TObjArray;
}
-//_____________________________________________________________________________
-
ClassImp(AliMUONhit)
fEloss=hits[7];
fPHfirst=(Int_t) hits[8];
fPHlast=(Int_t) hits[9];
-}
-ClassImp(AliMUONreccluster)
-
-ClassImp(AliMUONRecCluster)
-//_____________________________________________________________________
-AliMUONRecCluster::AliMUONRecCluster()
+ // modifs perso
+ fPTot=hits[10];
+ fCxHit=hits[11];
+ fCyHit=hits[12];
+ fCzHit=hits[13];
+}
+ClassImp(AliMUONcorrelation)
+//___________________________________________
+//_____________________________________________________________________________
+AliMUONcorrelation::AliMUONcorrelation(Int_t *idx, Float_t *x, Float_t *y)
{
- fDigits=0;
- fNdigit=-1;
+ //
+ // Creates a MUON correlation object
+ //
+ for(Int_t i=0; i<4; i++) {
+ fCorrelIndex[i] = idx[i];
+ fX[i] = x[i];
+ fY[i] = y[i];
+ }
+}
+ClassImp(AliMUONRawCluster)
+Int_t AliMUONRawCluster::Compare(TObject *obj)
+{
+ /*
+ AliMUONRawCluster *raw=(AliMUONRawCluster *)obj;
+ Float_t r=GetRadius();
+ Float_t ro=raw->GetRadius();
+ if (r>ro) return 1;
+ else if (r<ro) return -1;
+ else return 0;
+ */
+ AliMUONRawCluster *raw=(AliMUONRawCluster *)obj;
+ Float_t y=fY;
+ Float_t yo=raw->fY;
+ if (y>yo) return 1;
+ else if (y<yo) return -1;
+ else return 0;
+
}
-AliMUONRecCluster::
-AliMUONRecCluster(Int_t FirstDigit,Int_t Ichamber, Int_t Icathod)
+Int_t AliMUONRawCluster::
+BinarySearch(Float_t y, TArrayF coord, Int_t from, Int_t upto)
{
- fX = 0.;
- fY = 0.;
- fDigits = new TArrayI(10);
- fNdigit=0;
- AddDigit(FirstDigit);
- fChamber=Ichamber;
- fCathod=Icathod;
+ // Find object using a binary search. Array must first have been sorted.
+ // Search can be limited by setting upto to desired index.
+
+ Int_t low=from, high=upto-1, half;
+ while(high-low>1) {
+ half=(high+low)/2;
+ if(y>coord[half]) low=half;
+ else high=half;
+ }
+ return low;
}
-void AliMUONRecCluster::AddDigit(Int_t Digit)
+void AliMUONRawCluster::SortMin(Int_t *idx,Float_t *xdarray,Float_t *xarray,Float_t *yarray,Float_t *qarray, Int_t ntr)
{
- if (fNdigit==fDigits->GetSize()) {
- //enlarge the list by hand!
- Int_t *array= new Int_t[fNdigit*2];
- for(Int_t i=0;i<fNdigit;i++)
- array[i] = fDigits->At(i);
- fDigits->Adopt(fNdigit*2,array);
+ //
+ // Get the 3 closest points(cog) one can find on the second cathode
+ // starting from a given cog on first cathode
+ //
+
+ //
+ // Loop over deltax, only 3 times
+ //
+
+ Float_t xmin;
+ Int_t jmin;
+ Int_t id[3] = {-2,-2,-2};
+ Float_t jx[3] = {0.,0.,0.};
+ Float_t jy[3] = {0.,0.,0.};
+ Float_t jq[3] = {0.,0.,0.};
+ Int_t jid[3] = {-2,-2,-2};
+ Int_t i,j,imax;
+
+ if (ntr<3) imax=ntr;
+ else imax=3;
+ for(i=0;i<imax;i++){
+ xmin=1001.;
+ jmin=0;
+
+ for(j=0;j<ntr;j++){
+ if ((i == 1 && j == id[i-1])
+ ||(i == 2 && (j == id[i-1] || j == id[i-2]))) continue;
+ if (TMath::Abs(xdarray[j]) < xmin) {
+ xmin = TMath::Abs(xdarray[j]);
+ jmin=j;
+ }
+ } // j
+ if (xmin != 1001.) {
+ id[i]=jmin;
+ jx[i]=xarray[jmin];
+ jy[i]=yarray[jmin];
+ jq[i]=qarray[jmin];
+ jid[i]=idx[jmin];
+ }
+
+ } // i
+
+ for (i=0;i<3;i++){
+ if (jid[i] == -2) {
+ xarray[i]=1001.;
+ yarray[i]=1001.;
+ qarray[i]=1001.;
+ idx[i]=-1;
+ } else {
+ xarray[i]=jx[i];
+ yarray[i]=jy[i];
+ qarray[i]=jq[i];
+ idx[i]=jid[i];
+ }
}
- fDigits->AddAt(Digit,fNdigit);
- fNdigit++;
+
}
-AliMUONRecCluster::~AliMUONRecCluster()
+Int_t AliMUONRawCluster::PhysicsContribution()
{
- if (fDigits)
- delete fDigits;
+ Int_t iPhys=0;
+ Int_t iBg=0;
+ Int_t iMixed=0;
+ for (Int_t i=0; i<fMultiplicity; i++) {
+ if (fPhysicsMap[i]==2) iPhys++;
+ if (fPhysicsMap[i]==1) iMixed++;
+ if (fPhysicsMap[i]==0) iBg++;
+ }
+ if (iMixed==0 && iBg==0) {
+ return 2;
+ } else if ((iPhys != 0 && iBg !=0) || iMixed != 0) {
+ return 1;
+ } else {
+ return 0;
+ }
}
-Int_t AliMUONRecCluster::FirstDigitIndex()
-{
- fCurrentDigit=0;
- return fDigits->At(fCurrentDigit);
-}
+
+ClassImp(AliMUONreccluster)
+ClassImp(AliMUONsegmentation)
+ClassImp(AliMUONresponse)
+
+
+
+
+
-Int_t AliMUONRecCluster::NextDigitIndex()
-{
- fCurrentDigit++;
- if (fCurrentDigit<fNdigit)
- return fDigits->At(fCurrentDigit);
- else
- return InvalidDigitIndex();
-}
-Int_t AliMUONRecCluster::NDigits()
-{
- return fNdigit;
-}
-void AliMUONRecCluster::Finish()
-{
- // In order to reconstruct coordinates, one has to
- // get back to the digits which is not trivial here,
- // because we don't know where digits are stored!
- // Center Of Gravity, or other method should be
- // a property of AliMUON class!
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff