ITS/ITSRF.C

   1 #define __COMPILED__
   2 #ifdef __COMPILED__
   3 #include <iostream.h>
   4 #include <AliITS.h>
   5 #include <AliITSRiemannFit.h>
   6 #include <AliRun.h>
   7 #include <TROOT.h>
   8 #include <TClonesArray.h>
   9 #include <TFitter.h>
  10 #include <TMinuit.h>
  11 #include <TF2.h>
  12 #include <TFile.h>
  13 #include <TH1.h>
  14 #include <TH2.h>
  15 #include <TObjArray.h>
  16 #include <TParticle.h>
  17 #include <TTree.h>
  18 #include <TVector3.h>
  19 #include <TVirtualFitter.h>
  20 #endif
  21
  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)
  24 //
  25 void ITSRF(const char *filename="galice.root",
  26            const char *outfile="RiemannFit.dat",
  27            Int_t evnt=0) {
  28
  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);
  33
  34   // Get AliRun object from file or create it if not on file
  35   if (!gAlice) {
  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");
  39   }
  40
  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)
  45
  46   //
  47   //  Treat ITS information per module
  48   //
  49
  50   Int_t nmodules;
  51   ITS->InitModules(-1,nmodules);
  52   cout<<" Filling modules..."<<endl;
  53   ITS->FillModules(evnt,-1,nmodules," "," ");
  54
  55   //
  56   //  Get the recontruction tree
  57   //
  58   TTree *TR = gAlice->TreeR();
  59   Int_t nent = (Int_t)TR->GetEntries();
  60
  61   //
  62   // Initialize Riemann-fit class
  63   //
  64   AliITSRiemannFit *rfit = new AliITSRiemannFit();
  65   rfit->InitPoints(evnt,nent,ITS,TR,nparticles);
  66
  67   const TVector3 Bfield(0.0,0.0,0.2); // ALICE magnetic field
  68   //
  69   //   Setting the variables
  70   //
  71
  72   TParticle *MPart;
  73
  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;
  79
  80   FILE *outdat=fopen(outfile,"w");
  81
  82   cout<<endl<<" Looping on "<< nparticles <<" particles..."<<endl;
  83
  84   Int_t nAcc=0;
  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();
  89
  90     //
  91     //  get track fit
  92     //
  93     Pt = MPart->Pt();
  94     Px = MPart->Px();
  95     Py = MPart->Py();
  96     Pz = MPart->Pz();
  97     Ptot = MPart->P();
  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
 103     //
 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;
 112       }
 113       continue;
 114     }
 115
 116     //
 117     //
 118     //
 119     Int_t written=0;
 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);
 129     nAcc++;
 130
 131   } // loop on particles
 132   fclose(outdat);
 133
 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]);
 136
 137
 138   return;
 139 }
 140
 141