PWG1/TRD/macros/AliTrackletsinTRD.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ////////////////////////////////////////////////////////////////////////////
  17 //                                                                        //
  18 //  Reconstruction QA                                                     //
  19 //                                                                        //
  20 //  Authors: Yvonne C. Pachmayer <pachmay@physi.uni-heidelberg.de>        //
  21 //                                                                        //
  22 ////////////////////////////////////////////////////////////////////////////
  23
  24
  25 #define AliTrackletsinTRD_cxx
  26 #include "AliTrackletsinTRD.h"
  27 #include <TH2.h>
  28 #include <TStyle.h>
  29 #include <TCanvas.h>
  30
  31 void AliTrackletsinTRD::Loop()
  32 {
  33     //   In a ROOT session, you can do:
  34     //      Root > .L AliTrackletsinTRD.C
  35     //      Root > AliTrackletsinTRD t
  36     //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  37     //      Root > t.Show();       // Show values of entry 12
  38     //      Root > t.Show(16);     // Read and show values of entry 16
  39     //      Root > t.Loop();       // Loop on all entries
  40     //
  41
  42     //     This is the loop skeleton where:
  43     //    jentry is the global entry number in the chain
  44     //    ientry is the entry number in the current Tree
  45     //  Note that the argument to GetEntry must be:
  46     //    jentry for TChain::GetEntry
  47     //    ientry for TTree::GetEntry and TBranch::GetEntry
  48     //
  49     //       To read only selected branches, Insert statements like:
  50     // METHOD1:
  51     //    fChain->SetBranchStatus("*",0);  // disable all branches
  52     //    fChain->SetBranchStatus("branchname",1);  // activate branchname
  53     // METHOD2: replace line
  54     //    fChain->GetEntry(jentry);       //read all branches
  55     //by  b_branchname->GetEntry(ientry); //read only this branch
  56
  57     gROOT->SetStyle("Plain");
  58     gStyle->SetPalette(1);
  59
  60
  61     if (fChain == 0) return;
  62
  63     // open run loader and load gAlice, kinematics and header
  64     AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  65     if (!runLoader)
  66     {
  67         printf("Error: readKine:\nGetting run loader from file \"%s\" failed", filename);
  68         return;
  69     }
  70
  71     runLoader->LoadHeader();
  72     runLoader->LoadRecPoints("TRD");
  73     TObjArray *module = new TObjArray();
  74
  75     runLoader->CdGAFile();
  76
  77     Int_t nEvents = runLoader->GetNumberOfEvents();
  78     const Int_t nEventsarray=nEvents;
  79     AliTRDcluster *cls = 0;
  80
  81     Double_t x_clus[540][90];
  82     Double_t y_clus[540][90];
  83
  84     for(Int_t ev = 0; ev < nEvents - 1; ev++)
  85     {
  86         TTree *tree = runLoader->GetTreeR("TRD", 0);
  87         tree->SetBranchAddress("TRDcluster", &module);
  88
  89         int N = tree->GetEntries(); // number of chamber, max 540
  90         // Check number of clusters
  91         for(Int_t ind = 0; ind < N; ind++)
  92         {
  93             tree->GetEntry(ind);
  94             Int_t m = module->GetEntries();
  95
  96             for (Int_t j = 0; j < m; j++)    // loop over clusters of one chamber
  97             {
  98                 if(j>=90) break;
  99                 if (cls != 0) delete cls;
 100                 cls = (AliTRDcluster*)module->At(j);
 101                 x_clus[ind][j]=cls->GetX();
 102                 y_clus[ind][j]=cls->GetY();
 103
 104             }
 105         }
 106
 107
 108
 109         // loop over debug file and analysis
 110
 111         Float_t xbmin = 0., xbmax = 380.;
 112         /* // histogram for debugging purpose
 113         Int_t nxbins = 3000, nybins = 280;
 114         Float_t ybmin = -70., ybmax = 70.;
 115         hxy = new TH2F("hxy",";x;y",nxbins,xbmin,xbmax,nybins,ybmin,ybmax);
 116         */
 117
 118         Long64_t nentries = fChain->GetEntriesFast();
 119
 120         Long64_t nbytes = 0, nb = 0;
 121         Float_t xpos[6];
 122         Float_t ypos[6];
 123         Int_t counter_test;
 124         for (Long64_t jentry=0; jentry<nentries;jentry++) {
 125             Long64_t ientry = LoadTree(jentry);
 126             if (ientry < 0) break;
 127             nb = fChain->GetEntry(jentry);   nbytes += nb;
 128             xpos[flayer]=fxtrack;
 129             ypos[flayer]=fytrack;
 130
 131
 132             if(flayer==5)
 133             {
 134                 Float_t x[6]= {xpos[0],xpos[1],xpos[2],xpos[3],xpos[4],xpos[5]};
 135                 Float_t y[6]= {ypos[0],ypos[1],ypos[2],ypos[3],ypos[4],ypos[5]};
 136
 137                 gxy = new TGraph(6,x,y);
 138
 139                 TF1 *f1 = new TF1("f1","pol1",xbmin,xbmax);
 140                 gxy->Fit("f1","RQ");
 141                 // resulting function: y=ax+b --> x*f1->GetParameter(1)+f1->GetParameter(0)
 142
 143                 // graph and fit plotting only for checking
 144                 /*
 145                  c1 = new TCanvas("c1","hough Transform",0,0,800,800);
 146                  c1->cd(1);
 147                  hxy->Draw();
 148                  gxy->SetMarkerStyle(3);
 149                  gxy->Draw("p");
 150                  */
 151
 152                 Float_t px;
 153                 Float_t py;
 154
 155                 if(feventcounter==ev)
 156                 {
 157                     for(Int_t b=0;b<540;b++)
 158                     {
 159                         if(b==fdettracklet)
 160                         {
 161                             Int_t counter_distcalc;
 162                             Float_t distance_sum=0;
 163
 164                             for(Int_t c=0;c<90;c++)
 165                             {
 166                                 px=0;
 167                                 py=0;
 168
 169                                 px=x_clus[b][c];
 170                                 py=y_clus[b][c];
 171
 172                                 if(px!=0 &&  py!=0)
 173                                 {
 174                                     // Function to calculate the distance of a point to the above fitted function
 175                                     Double_t distance = 0;
 176                                     if((TMath::Sqrt(f1->GetParameter(1)*f1->GetParameter(1)+1))!=0) distance=TMath::Abs(f1->GetParameter(1)*px-py+f1->GetParameter(0))/(TMath::Sqrt(f1->GetParameter(1)*f1->GetParameter(1)+1));
 177                                     // if(distance<10) cout << eventcounter << " " << b << " " << c << " " <<   distance << endl;
 178                                     if(distance>0.6 && distance<=2)
 179                                     {
 180                                         counter_distcalc++;
 181                                         distance_sum+=distance;
 182                                     }
 183                                 }
 184
 185                             }
 186                             if(counter_distcalc!=0)
 187                             {
 188                                 // these are tracks which have additional clusters close by: 0.6 < distance <= 2
 189                             }
 190                             else
 191                             {
 192                                 // these are good tracks no additional clusters close by
 193                                cout <<" good tracks " <<  endl;
 194                             }
 195                             counter_distcalc=0;
 196                         }
 197                     }
 198                 }
 199
 200                 if(gxy) delete gxy;
 201                 if(f1) delete f1;
 202             } // end of 5 layer loop
 203
 204
 205
 206
 207
 208         }  // end loop over debug file and analysis
 209
 210         runLoader->GetNextEvent();
 211     }
 212
 213 }
 214