ITS/AlignITSmacro3Rphi.C

   1 //////////////////////////////////////////////////////////////////////////
   2 //  Alice ITS first detector alignment program.                         //
   3 //                                                                      //
   4 // version: 0.0.0 Draft.                                                //
   5 // Date: April 18 1999                                                  //
   6 // By: Bjorn S. Nilsen                                                  //
   7 //                                                                      //
   8 //////////////////////////////////////////////////////////////////////////
   9
  10 #include <fstream.h>
  11
  12 // Data structure to hold averaged clusts.
  13 struct ClustAl_sl{
  14     Int_t lay,lad,det;
  15     Float_t xg,yg,zg,xl,yl,zl;
  16 };
  17 struct ClustAl_tl{
  18     Int_t    track,nclust;  // track number and number of data points.
  19     ClustAl_sl *clust;      // data points to fit.
  20     Float_t  a,b,c,d,a0,b0,c0,d0,qual;  // fit parameters and fit quality.
  21     Float_t  px,py,pz,p,pt;
  22     // x=a+bz and y=c+dz;
  23     // x=a0+b0*y and z=c0+d0*y in coordinate system of clust[0].lay,lad,det
  24 };
  25
  26 //
  27 void AlignITSmacro3Rphi(const char *Rfilename="galice_ITS_B0.root",
  28                     const char *sfile="Align_ITS_B0",
  29                     Float_t sigma1=0.0,Float_t sigma2=0.0,Float_t sigma3=0.0,
  30                     Int_t evNumber=0) {
  31 /////////////////////////////////////////////////////////////////////////
  32 //   This macro is a small example of a ROOT macro illustrating how to
  33 //   read the output of GALICE and fill some histograms.
  34 //
  35 //     Root > .L AlignITSmacro2.C    // this loads the macro in memory
  36 //     Root > AlignITSmacro2();      // by default process first event
  37 //     Root > AlignITSmacro2("galice.root");   // process file galice.root
  38 //     Root > AlignITSmacro2("galice.root",3); // process file galice.root
  39 //                                                the first 4 events
  40 //       or
  41 //     Root > .x AlignITSmacro2.C;  //by default process first event
  42 //     Root > .x AlignITSmacro2.C("galice.root");   // process file galice.root
  43 //     Root > .x AlignITSmacro2.C("galice.root",3); // process file galice.root
  44 //                                                     the first 4 events
  45 /////////////////////////////////////////////////////////////////////////
  46 //
  47     gROOT->Reset();  // Reset root to it's default state
  48 // Dynamically link some shared libs
  49    if (gClassTable->GetID("AliRun") < 0) {
  50       gROOT->LoadMacro("loadlibs.C");
  51       loadlibs();
  52    } // end if gClassTable...
  53
  54    // Connect the Root Galice file containing Geometry, Kine and Clusts
  55    TFile *Rfile = (TFile*)gROOT->GetListOfFiles()->FindObject(Rfilename);
  56    if(!Rfile) Rfile = new TFile(Rfilename);
  57    printf("reading from file %s\n",Rfilename);
  58
  59    // Get AliRun object from file or create it if not on file
  60    if(!gAlice) {
  61       gAlice = (AliRun*)Rfile->Get("gAlice");
  62       if( gAlice) printf("AliRun object found on file\n");
  63       if(!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  64    } /* end if gAlice */
  65
  66    Float_t    v0[3] = {0.0,0.0,0.0};
  67    Float_t    tran[3] = {0.0,0.0,0.0},rot[3] = {0.0,0.0,0.0};
  68    Float_t    trans[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
  69                           2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
  70                           2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // cm
  71    Float_t    rots[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
  72                          2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
  73                          2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // rad
  74    Int_t      nparticles,Ntrkp,ntrk,Itimes,Isigmas;
  75    AliITS     *ITS = 0;
  76    TTree      *TH = 0;
  77    AliITSgeom *gm,gm2;
  78    char       Hfilename[80];
  79 //   char       Gfilename[80];
  80    TFile      *Hfile = 0;
  81    ClustAl_tl *trk = 0;
  82    Float_t    Rdta[6];
  83    Float_t    Fdta[24],Fdr[12],Fdta0[24],Fdr0[12];
  84    Int_t      Ndta[12];
  85    FILE       *fp1,*fp2;
  86
  87    for(Int_t evnt=0;evnt<=evNumber;evnt++){
  88       // define some variables for later use.
  89       nparticles = gAlice->GetEvent(evnt);
  90       printf("nparticles %d\n",nparticles);
  91       if (nparticles <= 0) continue; /* get next event */
  92
  93       // Get pointers to Alice detectors and Clusts containers
  94       ITS   = (AliITS*)gAlice->GetDetector("ITS");
  95       if(!ITS) return;          /* error no ITS data exit */
  96       TH    = gAlice->TreeH();
  97       Ntrkp = TH->GetEntries();
  98       gm    = ITS->GetITSgeom();
  99
 100       // Array (stucture) of clusts for the first and second layer
 101       // this should be replaced with either clusters or digits
 102       // when they are proporly defined.
 103       trk = new ClustAl_tl[Ntrkp];
 104
 105       printf("Ntrkp=%d\n",Ntrkp);
 106
 107       HitsToClustAl(trk,ntrk,Ntrkp,TH,ITS,1.0);
 108       printf("Filled data structures ntrk=%d fitting lines next\n",ntrk);
 109 //
 110       for(Itimes=1;Itimes<2;Itimes++){
 111           if(Itimes==0){
 112               fp1 = fopen("RvariationsR.csv","w");
 113               fp2 = fopen("RvariationsN.csv","w");
 114             fprintf(fp1,"Rr,Rs,Rrx1,Rerx1,Rsrx1,Rserx1,Rrz1,Rerz1,Rsrz1,Rserz1,"
 115                     "Rrd1,Rerd1,Rsrd1,Resrd1,"
 116                              "Rrx2,Rerx2,Rsrx2,Rserx2,Rrz2,Rerz2,Rsrz2,Rserz2,"
 117                     "Rrd2,Rerd2,Rsrd2,Resrd2,"
 118                              "Rrx3,Rerx3,Rsrx3,Rserx3,Rrz3,Rerz3,Rsrz3,Rserz3,"
 119                     "Rrd3,Rerd3,Rsrd3,Resrd3,"
 120                              "Rrx4,Rerx4,Rsrx4,Rserx4,Rrz4,Rerz4,Rsrz4,Rserz4,"
 121                     "Rrd4,Rerd4,Rsrd4,Resrd4,"
 122                              "Rrx5,Rerx5,Rsrx5,Rserx5,Rrz5,Rerz5,Rsrz5,Rserz5,"
 123                     "Rrd5,Rerd5,Rsrd5,Resrd5,"
 124                             "Rrx6,Rerx6,Rsrx6,Rserx6,Rrz6,Rerz6,Rsrz6,Rserz6,"
 125                     "Rrd6,Rerd6,Rsrd6,Resrd6");
 126             fprintf(fp2,"Rr,Rs,RN1,RN2,RN3,RN4,RN5,RN6,RN7,RN8,RN9,RN10,"
 127                     "RN50,Rcut");
 128           }else if(Itimes==1){
 129               fp1 = fopen("RPhivariationsR.csv","w");
 130               fp2 = fopen("RPhivariationsN.csv","w");
 131             fprintf(fp1,"Fr,Fs,Frx1,Ferx1,Fsrx1,Fserx1,Frz1,Ferz1,Fsrz1,Fserz1,"
 132                     "Frd1,Ferd1,Fsrd1,Fesrd1,"
 133                              "Frx2,Ferx2,Fsrx2,Fserx2,Frz2,Ferz2,Fsrz2,Fserz2,"
 134                     "Frd2,Ferd2,Fsrd2,Fesrd2,"
 135                              "Frx3,Ferx3,Fsrx3,Fserx3,Frz3,Ferz3,Fsrz3,Fserz3,"
 136                     "Frd3,Ferd3,Fsrd3,Fesrd3,"
 137                              "Frx4,Ferx4,Fsrx4,Fserx4,Frz4,Ferz4,Fsrz4,Fserz4,"
 138                     "Frd4,Ferd4,Fsrd4,Fesrd4,"
 139                              "Frx5,Ferx5,Fsrx5,Fserx5,Frz5,Ferz5,Fsrz5,Fserz5,"
 140                     "Frd5,Ferd5,Fsrd5,Fesrd5,"
 141                             "Frx6,Ferx6,Fsrx6,Fserx6,Frz6,Ferz6,Fsrz6,Fserz6,"
 142                     "Frd6,Ferd6,Fsrd6,Fesrd6");
 143             fprintf(fp2,"Fr,Fs,FN1,FN2,FN3,FN4,FN5,FN6,FN7,FN8,FN9,FN10,"
 144                     "FN50,Fcut");
 145           }else if(Itimes==2){
 146               fp1 = fopen("ZvariationsR.csv","w");
 147               fp2 = fopen("ZvariationsN.csv","w");
 148             fprintf(fp1,"Zr,Zs,Zrx1,Zerx1,Zsrx1,Zserx1,Zrz1,Zerz1,Zsrz1,Zserz1,"
 149                     "Zrd1,Zerd1,Zsrd1,Zesrd1,"
 150                              "Zrx2,Zerx2,Zsrx2,Zserx2,Zrz2,Zerz2,Zsrz2,Zserz2,"
 151                     "Zrd2,Zerd2,Zsrd2,Zesrd2,"
 152                              "Zrx3,Zerx3,Zsrx3,Zserx3,Zrz3,Zerz3,Zsrz3,Zserz3,"
 153                     "Zrd3,Zerd3,Zsrd3,Zesrd3,"
 154                              "Zrx4,Zerx4,Zsrx4,Zserx4,Zrz4,Zerz4,Zsrz4,Zserz4,"
 155                     "Zrd4,Zerd4,Zsrd4,Zesrd4,"
 156                              "Zrx5,Zerx5,Zsrx5,Zserx5,Zrz5,Zerz5,Zsrz5,Zserz5,"
 157                     "Zrd5,Zerd5,Zsrd5,Zesrd5,"
 158                             "Zrx6,Zerx6,Zsrx6,Zserx6,Zrz6,Zerz6,Zsrz6,Zserz6,"
 159                     "Zrd6,Zerd6,Zsrd6,Zesrd6");
 160             fprintf(fp2,"Zr,Zs,ZN1,ZN2,ZN3,ZN4,ZN5,ZN6,ZN7,ZN8,ZN9,ZN10,"
 161                     "ZN50,Zcut");
 162           }else if(Itimes==3){
 163               fp1 = fopen("AvariationsR.csv","w");
 164               fp2 = fopen("AvariationsN.csv","w");
 165             fprintf(fp1,"Ar,As,Arx1,Aerx1,Asrx1,Aserx1,Arz1,Aerz1,Asrz1,Aserz1,"
 166                     "Ard1,Aerd1,Asrd1,Aesrd1,"
 167                              "Arx2,Aerx2,Asrx2,Aserx2,Arz2,Aerz2,Asrz2,Aserz2,"
 168                     "Ard2,Aerd2,Asrd2,Aesrd2,"
 169                              "Arx3,Aerx3,Asrx3,Aserx3,Arz3,Aerz3,Asrz3,Aserz3,"
 170                     "Ard3,Aerd3,Asrd3,Aesrd3,"
 171                              "Arx4,Aerx4,Asrx4,Aserx4,Arz4,Aerz4,Asrz4,Aserz4,"
 172                     "Ard4,Aerd4,Asrd4,Aesrd4,"
 173                              "Arx5,Aerx5,Asrx5,Aserx5,Arz5,Aerz5,Asrz5,Aserz5,"
 174                     "Ard5,Aerd5,Asrd5,Aesrd5,"
 175                             "Arx6,Aerx6,Asrx6,Aserx6,Arz6,Aerz6,Asrz6,Aserz6,"
 176                     "Ard6,Aerd6,Asrd6,Aesrd6");
 177             fprintf(fp2,"Ar,As,AN1,AN2,AN3,AN4,AN5,AN6,AN7,AN8,AN9,AN10,"
 178                     "AN50,Acut");
 179           }else if(Itimes==4){
 180               fp1 = fopen("BvariationsR.csv","w");
 181               fp2 = fopen("BvariationsN.csv","w");
 182             fprintf(fp1,"Br,Bs,Brx1,Berx1,Bsrx1,Bserx1,Brz1,Berz1,Bsrz1,Bserz1,"
 183                     "Brd1,Berd1,Bsrd1,Besrd1,"
 184                              "Brx2,Berx2,Bsrx2,Bserx2,Brz2,Berz2,Bsrz2,Bserz2,"
 185                     "Brd2,Berd2,Bsrd2,Besrd2,"
 186                              "Brx3,Berx3,Bsrx3,Bserx3,Brz3,Berz3,Bsrz3,Bserz3,"
 187                     "Brd3,Berd3,Bsrd3,Besrd3,"
 188                              "Brx4,Berx4,Bsrx4,Bserx4,Brz4,Berz4,Bsrz4,Bserz4,"
 189                     "Brd4,Berd4,Bsrd4,Besrd4,"
 190                              "Brx5,Berx5,Bsrx5,Bserx5,Brz5,Berz5,Bsrz5,Bserz5,"
 191                     "Brd5,Berd5,Bsrd5,Besrd5,"
 192                             "Brx6,Berx6,Bsrx6,Bserx6,Brz6,Berz6,Bsrz6,Bserz6,"
 193                     "Brd6,Berd6,Bsrd6,Besrd6");
 194             fprintf(fp2,"Br,Bs,BN1,BN2,BN3,BN4,BN5,BN6,BN7,BN8,BN9,BN10,"
 195                     "BN50,Bcut");
 196           }else if(Itimes==5){
 197               fp1 = fopen("CvariationsR.csv","w");
 198               fp2 = fopen("CvariationsN.csv","w");
 199             fprintf(fp1,"Cr,Cs,Crx1,Cerx1,Csrx1,Cserx1,Crz1,Cerz1,Csrz1,Cserz1,"
 200                     "Crd1,Cerd1,Csrd1,Cesrd1,"
 201                              "Crx2,Cerx2,Csrx2,Cserx2,Crz2,Cerz2,Csrz2,Cserz2,"
 202                     "Crd2,Cerd2,Csrd2,Cesrd2,"
 203                              "Crx3,Cerx3,Csrx3,Cserx3,Crz3,Cerz3,Csrz3,Cserz3,"
 204                     "Crd3,Cerd3,Csrd3,Cesrd3,"
 205                              "Crx4,Cerx4,Csrx4,Cserx4,Crz4,Cerz4,Csrz4,Cserz4,"
 206                     "Crd4,Cerd4,Csrd4,Cesrd4,"
 207                              "Crx5,Cerx5,Csrx5,Cserx5,Crz5,Cerz5,Csrz5,Cserz5,"
 208                     "Crd5,Cerd5,Csrd5,Cesrd5,"
 209                     "Crx6,Cerx6,Csrx6,Cserx6,Crz6,Cerz6,Csrz6,Cserz6,"
 210                     "Crd6,Cerd6,Csrd6,Cesrd6");
 211             fprintf(fp2,"Cr,Cs,CN1,CN2,CN3,CN4,CN5,CN6,CN7,CN8,CN9,CN10,"
 212                     "CN50,Ccut");
 213           } // end if
 214           fprintf(fp1,"\n");
 215           fprintf(fp2,"\n");
 216           for(Isigmas=0;Isigmas<15;Isigmas++){
 217 //
 218 //        tran[0] = sigma1;
 219 //        tran[1] = sigma2;
 220 //        tran[2] = sigma3;
 221           if(Itimes==0){ tran[0] = trans[Isigmas];
 222           }else tran[0] = 0.0;
 223           if(Itimes==1){ tran[1] = trans[Isigmas];
 224           }else tran[1] = 0.0;
 225           if(Itimes==2){ tran[2] = trans[Isigmas];
 226           }else tran[2] = 0.0;
 227           if(Itimes==3){ rot[0] = rots[Isigmas];
 228           }else rot[0] = 0.0;
 229           if(Itimes==4){ rot[1] = rots[Isigmas];
 230           }else rot[1] = 0.0;
 231           if(Itimes==5){ rot[2] = rots[Isigmas];
 232           }else rot[2] = 0.0;
 233           printf("tran= %e %e %e (cm), rot=%e %e %e (rad)\n",
 234                  tran[0],tran[1],tran[2],rot[0],rot[1],rot[2]);
 235 //
 236           gm2 = *gm;
 237           gm2.RandomCylindericalChange(tran,rot);
 238
 239           FillGlobalPositions(trk,ntrk,&gm2);
 240           // setup to save all created histograms.
 241           sprintf(Hfilename,"%s_%04.0fr%04.0fp%04.0fz%04.0fx%04.0fy%04.0fz.root",
 242                   sfile,
 243                   10000.*tran[0],10000.*tran[1],10000.*tran[2],
 244                   10000.* rot[0],10000.* rot[1],10000.* rot[2]);
 245           Hfile = (TFile*)gROOT->GetListOfFiles()->FindObject(Hfilename);
 246           if(Hfile) Hfile->Close();
 247           Hfile = new TFile(Hfilename,"RECREATE","Histograms from AlignITS.C");
 248           printf("Histograms saved to file %s\n",Hfilename);
 249           //
 250           PlotGeomChanges(gm,&gm2,Hfile,Rdta);
 251           //
 252           // fit all tracks and do a track quality hist.
 253           FitAllTracks(trk,ntrk,v0,&gm2,sfile,Hfile,Fdta,Ndta);
 254 //
 255           for(Int_t l=0;l<6;l++){
 256               Fdr[2*l+0] = sqrt(Fdta[4*l+0]*Fdta[4*l+0] +
 257                                 Fdta[4*l+2]*Fdta[4*l+2]);
 258               Fdr[2*l+1] = sqrt(Fdta[4*l+0]*Fdta[4*l+0]*
 259                                 Fdta[4*l+1]*Fdta[4*l+1] +
 260                                 Fdta[4*l+2]*Fdta[4*l+2]*
 261                                 Fdta[4*l+3]*Fdta[4*l+3])/Fdr[2*l+0];
 262           } // end for l
 263           if(Isigmas==0){
 264               for(Int_t fp=0;fp<24;fp++){
 265                   Fdta0[fp] = Fdta[fp];
 266                   if(fp<12) Fdr0[fp]  = Fdr[fp];
 267               } // end for fp
 268           } // end if Itimes==0&&Isigmas==0
 269           if(Itimes==0) fprintf(fp1,"%e,%e,",tran[0],Rdta[0]);
 270           else if(Itimes==1) fprintf(fp1,"%e,%e,",tran[1],Rdta[1]);
 271           else if(Itimes==2) fprintf(fp1,"%e,%e,",tran[2],Rdta[2]);
 272           else if(Itimes==3) fprintf(fp1,"%e,%e,",rot[0],Rdta[3]);
 273           else if(Itimes==4) fprintf(fp1,"%e,%e,",rot[1],Rdta[4]);
 274           else if(Itimes==5) fprintf(fp1,"%e,%e,",rot[2],Rdta[5]);
 275           if(Itimes==0) fprintf(fp2,"%e,%e,",tran[0],Rdta[0]);
 276           else if(Itimes==1) fprintf(fp2,"%e,%e,",tran[1],Rdta[1]);
 277           else if(Itimes==2) fprintf(fp2,"%e,%e,",tran[2],Rdta[2]);
 278           else if(Itimes==3) fprintf(fp2,"%e,%e,",rot[0],Rdta[3]);
 279           else if(Itimes==4) fprintf(fp2,"%e,%e,",rot[1],Rdta[4]);
 280           else if(Itimes==5) fprintf(fp2,"%e,%e,",rot[2],Rdta[5]);
 281           fprintf(fp1,"%e,%e,%e,%e,",Fdta[0],Fdta[1],Fdta[0]/Fdta0[0],Fdta[1]/Fdta0[0]);
 282           fprintf(fp1,"%e,%e,%e,%e,",Fdta[2],Fdta[3],Fdta[2]/Fdta0[2],Fdta[3]/Fdta0[2]);
 283           fprintf(fp1,"%e,%e,%e,%e,",Fdr[0],Fdr[1],Fdr[0]/Fdr0[0],Fdr[1]/Fdr0[0]);
 284           fprintf(fp1,"%e,%e,%e,%e,",Fdta[4],Fdta[5],Fdta[4]/Fdta0[4],Fdta[5]/Fdta0[4]);
 285           fprintf(fp1,"%e,%e,%e,%e,",Fdta[6],Fdta[7],Fdta[6]/Fdta0[6],Fdta[7]/Fdta0[6]);
 286           fprintf(fp1,"%e,%e,%e,%e,",Fdr[2],Fdr[3],Fdr[2]/Fdr0[2],Fdr[3]/Fdr0[2]);
 287           fprintf(fp1,"%e,%e,%e,%e,",Fdta[8],Fdta[9],Fdta[8]/Fdta0[8],Fdta[9]/Fdta0[8]);
 288           fprintf(fp1,"%e,%e,%e,%e,",Fdta[10],Fdta[11],Fdta[10]/Fdta0[10],Fdta[11]/Fdta0[10]);
 289           fprintf(fp1,"%e,%e,%e,%e,",Fdr[4],Fdr[5],Fdr[4]/Fdr0[4],Fdr[5]/Fdr0[4]);
 290           fprintf(fp1,"%e,%e,%e,%e,",Fdta[12],Fdta[12],Fdta[12]/Fdta0[12],Fdta[13]/Fdta0[12]);
 291           fprintf(fp1,"%e,%e,%e,%e,",Fdta[14],Fdta[13],Fdta[14]/Fdta0[14],Fdta[15]/Fdta0[14]);
 292           fprintf(fp1,"%e,%e,%e,%e,",Fdr[6],Fdr[7],Fdr[6]/Fdr0[6],Fdr[7]/Fdr0[6]);
 293           fprintf(fp1,"%e,%e,%e,%e,",Fdta[16],Fdta[15],Fdta[16]/Fdta0[16],Fdta[17]/Fdta0[16]);
 294           fprintf(fp1,"%e,%e,%e,%e,",Fdta[18],Fdta[17],Fdta[18]/Fdta0[18],Fdta[19]/Fdta0[18]);
 295           fprintf(fp1,"%e,%e,%e,%e,",Fdr[8],Fdr[9],Fdr[8]/Fdr0[8],Fdr[9]/Fdr0[8]);
 296           fprintf(fp1,"%e,%e,%e,%e,",Fdta[20],Fdta[19],Fdta[20]/Fdta0[20],Fdta[21]/Fdta0[20]);
 297           fprintf(fp1,"%e,%e,%e,%e,",Fdta[22],Fdta[21],Fdta[22]/Fdta0[22],Fdta[23]/Fdta0[22]);
 298           fprintf(fp1,"%e,%e,%e,%e\n",Fdr[10],Fdr[11],Fdr[10]/Fdr0[10],Fdr[11]/Fdr0[10]);
 299           fprintf(fp2,"%d,%d,%d,%d,",Ndta[0],Ndta[1],Ndta[2],Ndta[3]);
 300           fprintf(fp2,"%d,%d,%d,%d,",Ndta[4],Ndta[5],Ndta[6],Ndta[7]);
 301           fprintf(fp2,"%d,%d,%d,%d\n",Ndta[8],Ndta[9],Ndta[10],Ndta[11]);
 302 //        FitVertexAll(trk,ntrk,sfile,Hfile);
 303           // Find all 2 track vertecies and hist. them.
 304           Hfile->Close();
 305           //
 306       } // end for Isigmas
 307           fclose(fp1);
 308           fclose(fp2);
 309       } // end for Itimes
 310 //
 311       printf("Event %d done\n",evnt);
 312 //
 313       deleteClustAl(trk,ntrk); // subrotine to delet memory allocated
 314                                // inside HitsToclustAl.
 315       delete[] trk;            // now delet memory allocated above.
 316    } // end for evnt
 317    Rfile->Close();
 318    return;
 319 }