5 #include <AliITSRiemannFit.h>
8 #include <TClonesArray.h>
15 #include <TObjArray.h>
16 #include <TParticle.h>
19 #include <TVirtualFitter.h>
22 // Writes the data of helices and the position of the true II-ary, to perform fit
23 // of II-ary and calculation of resolution (in another macro)
25 void ITSRF(const char *filename="galice.root",
26 const char *outfile="RiemannFit.dat",
29 // Connect the Root Galice file containing Geometry, Kine and Hits
30 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
31 if (!file){ file = new TFile(filename);printf("Opening new file\n");}
32 printf("Root input file is %s\n",filename);
34 // Get AliRun object from file or create it if not on file
36 gAlice = (AliRun*)file->Get("gAlice");
37 if (gAlice) printf("AliRun object found on file\n");
38 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
41 // Get pointers to Alice detectors and Hits containers
42 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
43 if(!ITS) {cout<<"not ITS"<<endl; return;}
44 Int_t nparticles = gAlice->GetEvent(evnt); // number of particles generated (including daughters)
47 // Treat ITS information per module
51 ITS->InitModules(-1,nmodules);
52 cout<<" Filling modules..."<<endl;
53 ITS->FillModules(evnt,-1,nmodules," "," ");
56 // Get the recontruction tree
58 TTree *TR = gAlice->TreeR();
59 Int_t nent = (Int_t)TR->GetEntries();
62 // Initialize Riemann-fit class
64 AliITSRiemannFit *rfit = new AliITSRiemannFit();
65 rfit->InitPoints(evnt,nent,ITS,TR,nparticles);
67 const TVector3 Bfield(0.0,0.0,0.2); // ALICE magnetic field
69 // Setting the variables
74 Double_t r0[3]={0.0,0.0,0.0}, w, chi;
75 Double_t Pt, Px, Py, Pz, charge, Ptot;
76 Double_t x0, y0, rho, fd0, fphi, z0, vpar, chisql;
77 Double_t CorrLin, ResSum;
78 Int_t status[4]={0,0,0,0}, gfit_status, PdgPart;
80 FILE *outdat=fopen(outfile,"w");
82 cout<<endl<<" Looping on "<< nparticles <<" particles..."<<endl;
85 for(Int_t part = 0; part < nparticles ; part++) {
86 if (!(part%100)) cout<<part<<"\r";
87 MPart = (TParticle*)gAlice->Particle(part);
88 PdgPart = MPart->GetPdgCode();
98 charge = PdgPart/TMath::Abs(PdgPart); // +/- 1
99 gfit_status = (Int_t)rfit->FitHelix(part, charge, Px, Py, Pz, fd0, fphi, x0, y0,
100 rho, w, z0, vpar, chisql, CorrLin, ResSum);
101 // With this offset chi = 0 at the
102 // distance of closest approach in bending plane
104 chi=-TMath::Pi()-fphi;
105 if(gfit_status > 0) {
106 switch(gfit_status) {
107 case 1: status[0]++;break;
108 case 2: status[1]++;break;
109 case 11: status[2]++;break;
110 case 12: status[3]++; break;
111 default: cout<<"CASO NON PREVISTO IN K: "<<gfit_status<<endl;break;
120 written=fprintf(outdat,"%d ", part);
121 written=fprintf(outdat,"%f ", r0[0]);
122 written=fprintf(outdat,"%f ", r0[1]);
123 written=fprintf(outdat,"%f \n", r0[2]);
124 written=fprintf(outdat,"%f \n", ResSum);
125 written=fprintf(outdat," %f ", rho);
126 written=fprintf(outdat,"%f ", w);
127 written=fprintf(outdat,"%f ", fphi);
128 written=fprintf(outdat,"%f \n", vpar);
131 } // loop on particles
134 printf(" missed fits because of (resp.):\n no six points, no solution to cubic equation, \n z-direction inversion, tangent limit \n");
135 printf(" %d %d %d %d %d \n",status[0],status[1],status[2],status[3]);