TRD/AliTRDclusterErrors.C

   1 #ifndef __CINT__
   2   #include <iostream.h>
   3
   4   #include "AliTRDtracker.h"
   5   #include "AliTRDcluster.h"
   6   #include "AliTRDhit.h"
   7   #include "AliTRDv1.h"
   8   #include "AliTRDgeometry.h"
   9   #include "AliTRDparameter.h"
  10
  11   #include "alles.h"
  12   #include "TFile.h"
  13   #include "TStopwatch.h"
  14
  15 #endif
  16
  17 void AliTRDclusterErrors() {
  18
  19   Int_t No_of_tracks_to_analyze = 40;
  20
  21   const Int_t tbpp = 15;
  22   const Int_t nPlanes = 6;
  23   const Int_t ntb = tbpp * nPlanes;
  24
  25   TH1F *hy = new TH1F("delta R*phi","Cluster displacement in R*phi",200,-1.,1.);
  26   TH1F *hyp = new TH1F("delta R*phi pos","delta R*phi, positive",200,-1.,1.);
  27   TH1F *hym = new TH1F("delta R*phi neg","delta R*phi, negative",200,-1.,1.);
  28
  29   TH1F *hyn = new TH1F("Norm., d(R*phi)","Norm. cluster displacement in R*phi",400,-8.,8.);
  30   TH1F *hz = new TH1F("delta Z","Cluster displacement in Z",300,-10.,50.);
  31   TH2F *hy2 = new TH2F("Amp vs delta R*phi","Amplitude versus delta R*phi",200,-5.,5.,200,0.,600.);
  32   TH2F *herr = new TH2F("sigmaY vs delta R*phi","sigmaY vs delta R*phi",200,-1,1,200,0.,0.1);
  33   TH2F *hy3 = new TH2F("Position within pad vs delta R*phi","Position within pad vs delta R*phi",200,-1.,1.,200,-0.5,1.5);
  34   TH2F *hy4 = new TH2F("local tb vs delta R*phi","local tb vs delta R*phi",200,-1.,1.,20,-2.5,17.5);
  35
  36   hy->SetXTitle("Displacement, cm");
  37   hyn->SetXTitle("Displacement, SigmaY");
  38   hy2->SetXTitle("Displacement, cm");
  39   hy2->SetYTitle("Amplitude");
  40   hz->SetXTitle("Displacement, cm");
  41   hy3->SetXTitle("Displacement, cm");
  42   hy3->SetYTitle("Position, cm");
  43   hy4->SetXTitle("Displacement, cm");
  44   hy4->SetYTitle("local time bin");
  45
  46   /*
  47   // Dynamically link some shared libs
  48   if (gClassTable->GetID("AliRun") < 0) {
  49     gROOT->LoadMacro("loadlibs.C");
  50     loadlibs();
  51     cout << "Loaded shared libraries" << endl;
  52   }
  53   */
  54
  55   // Load clusters
  56   Char_t *clusterfile = "AliTRDclusters.root";
  57   Int_t   nEvent  = 0;
  58
  59   TObjArray carray(2000);
  60   TObjArray *ClustersArray = &carray;
  61   TFile *geofile =TFile::Open("AliTRDclusters.root");
  62
  63   AliTRDparameter *par = (AliTRDparameter*) geofile->Get("TRDparameter");
  64   AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry");
  65
  66   AliTRDtracker *Tracker = new AliTRDtracker(geofile);
  67   Tracker->SetEventNumber(nEvent);
  68   Tracker->ReadClusters(ClustersArray,clusterfile);
  69   Int_t nClusters = carray.GetEntriesFast();
  70
  71   printf("Total number of clusters %d \n", nClusters);
  72
  73   Char_t *alifile = "galice.root";
  74   // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
  75   TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
  76   if (!gafl) {
  77     cout << "Open the ALIROOT-file " << alifile << endl;
  78     gafl = new TFile(alifile);
  79   }
  80   else {
  81     cout << alifile << " is already open" << endl;
  82   }
  83
  84   // Get AliRun object from file or create it if not on file
  85   gAlice = (AliRun*) gafl->Get("gAlice");
  86   if (gAlice)
  87     cout << "AliRun object found on file" << endl;
  88   else
  89     gAlice = new AliRun("gAlice","Alice test program");
  90
  91   AliTRDv1       *fTRD           = (AliTRDv1*) gAlice->GetDetector("TRD");
  92
  93   // Import the Trees for the event nEvent in the file
  94   Int_t nparticles = gAlice->GetEvent(nEvent);
  95
  96   TObjArray *particles=gAlice->Particles();
  97
  98   TTree *hitTree = gAlice->TreeH();
  99   Int_t  nBytes = 0;
 100
 101   // Get the number of entries in the hit tree
 102   // (Number of primary particles creating a hit somewhere)
 103   Int_t nTrack = (Int_t) hitTree->GetEntries();
 104   cout << " Found " << nTrack << " primary particles with hits" << endl;
 105   No_of_tracks_to_analyze = TMath::Min(No_of_tracks_to_analyze, nTrack);
 106
 107
 108   // Loop through particles and fill histoes
 109
 110   Float_t hitY[ntb];
 111   Float_t hitZ[ntb];
 112   Float_t hitO[ntb];
 113   Float_t clusterY[ntb];
 114   Float_t clusterZ[ntb];
 115   Float_t clusterQ[ntb];
 116   Float_t clusterSigmaY[ntb];
 117   Float_t pos[3];
 118   Float_t rot[3];
 119   Float_t global[3];
 120   Int_t det = 0;
 121   Int_t track_index[3];
 122
 123   printf("\n");
 124
 125   for (Int_t ii=0; ii<No_of_tracks_to_analyze; ii++) {
 126     printf("track %d out of %d \n", ii+1 , No_of_tracks_to_analyze);
 127
 128     for(Int_t plane = 0; plane < nPlanes; plane++) {
 129       for(Int_t ltb = 14; ltb > -1; ltb--) {
 130
 131         if(ii >= nTrack) continue;
 132
 133         TParticle *p = gAlice->Particle(ii);
 134         if (p->GetFirstMother()>=0) continue;
 135         TParticlePDG *pdg = p->GetPDG();
 136         Float_t charge=pdg->Charge();
 137         if(TMath::Abs(charge) < 0.5) continue;
 138
 139         Int_t gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
 140         Double_t x = Tracker->GetX(0,plane,ltb);
 141
 142         // loop through clusters
 143
 144         Bool_t cluster_found = kFALSE;
 145         Int_t nw = 0;
 146
 147         for (Int_t i = 0; i < nClusters; i++) {
 148
 149           AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);
 150
 151           nw = cl->GetDetector();
 152
 153           nw = fGeo->GetPlane(nw);
 154
 155           if(nw != plane) continue;
 156           for(Int_t j=0; j<3; j++) track_index[j] = cl->GetLabel(j);
 157           if((track_index[0] != ii) &&
 158              (track_index[1] != ii) &&
 159              (track_index[2] != ii)) continue;
 160
 161           nw = cl->GetLocalTimeBin();
 162           if(nw != ltb) continue;
 163
 164           clusterY[gtb] = cl->GetY();
 165           clusterZ[gtb] = cl->GetZ();
 166           clusterQ[gtb] = cl->GetQ();
 167
 168           clusterSigmaY[gtb] = TMath::Sqrt(cl->GetSigmaY2());
 169           cluster_found = kTRUE;
 170           break;
 171         }
 172
 173         if(!cluster_found) continue;
 174
 175         gAlice->ResetHits();
 176
 177         //      nBytes += hitTree->GetEvent(nPrimaries - ii - 1);
 178         nBytes += hitTree->GetEvent(nTrack - ii - 1);
 179
 180         // Loop through the TRD hits
 181         Bool_t found_hit = kFALSE;
 182
 183
 184         for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
 185             hit;
 186             hit = (AliTRDhit *) fTRD->NextHit()) {
 187           nw = hit->Track();
 188           if(nw != ii) continue;
 189           det   = hit->GetDetector();
 190           nw = fGeo->GetPlane(det);
 191           if(nw != plane) continue;
 192
 193
 194           pos[0]=hit->X();
 195           pos[1]=hit->Y();
 196           pos[2]=hit->Z();
 197           fGeo->Rotate(det,pos,rot);
 198
 199           if(TMath::Abs(rot[0]-x) > 0.01) continue;
 200           hitY[gtb] = rot[1];
 201           hitZ[gtb] = rot[2];
 202
 203           Float_t col0        = par->GetCol0(plane);
 204           Float_t colPadSize  = par->GetColPadSize(plane);
 205           Float_t colH = (Int_t ((rot[1] -  col0)/colPadSize)) * colPadSize;
 206           hitO[gtb] = (rot[1] -  col0) - colH;
 207           found_hit = kTRUE;
 208           break;
 209         }
 210
 211         if(!found_hit) continue;
 212
 213         /*
 214         printf("gtb: %d, x: %f, rot[0]: %f, Yhit: %f, Ycl: %f\n",
 215                gtb, x, rot[0], rot[1], clusterY[gtb]);
 216         printf("\n                            Zhit - Zcl = %f - %f = %f\n",
 217                rot[2], clusterZ[gtb], rot[2] - clusterZ[gtb]);
 218
 219
 220
 221         printf("found hit within dx = %f - %f \n",rot[0],x);
 222         printf("pos: %f, %f, %f \n",pos[0],pos[1],pos[2]);
 223         printf("rot: %f, %f, %f \n",rot[0],rot[1],rot[2]);
 224         printf("cluster: %d, %f, %f \n",gtb,clusterY[gtb],clusterZ[gtb]);
 225         */
 226
 227
 228         hy->Fill(hitY[gtb]-clusterY[gtb]);
 229         if(charge > 0) hyp->Fill(hitY[gtb]-clusterY[gtb]);
 230         else if(charge < 0) hym->Fill(hitY[gtb]-clusterY[gtb]);
 231
 232         if((clusterQ[gtb]>10)&&(clusterSigmaY[gtb]>0))
 233           hyn->Fill((hitY[gtb]-clusterY[gtb])/clusterSigmaY[gtb]);
 234         hz->Fill(hitZ[gtb]-clusterZ[gtb]);
 235         hy2->Fill(hitY[gtb]-clusterY[gtb],clusterQ[gtb]);
 236         hy3->Fill(hitY[gtb]-clusterY[gtb],hitO[gtb]);
 237         hy4->Fill(hitY[gtb]-clusterY[gtb],(Float_t)(tbpp - 1 - gtb%tbpp));
 238         herr->Fill(hitY[gtb]-clusterY[gtb],clusterSigmaY[gtb]);
 239       }
 240     }
 241   }
 242
 243   gStyle->SetOptStat(1);
 244   gStyle->SetOptFit(1);
 245
 246   TCanvas* c = new TCanvas("c", "c", 110, 110, 810, 840);
 247   c->SetFillColor(10);
 248   c->Divide(2,2);
 249   c->cd(1); hy->SetLineWidth(2); hy->SetFillColor(29); hy->Draw();
 250   c->cd(2); hz->SetLineWidth(2); hz->SetFillColor(29); hz->Draw();
 251   c->cd(3); hyn->Draw();
 252   c->cd(4); hy4->Draw();
 253
 254   TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940);
 255   c1->SetFillColor(10);
 256   c1->Divide(2,2);
 257   c1->cd(1); hyp->SetLineWidth(2); hyp->SetFillColor(29); hyp->Fit("gaus");
 258   c1->cd(2); hym->SetLineWidth(2); hym->SetFillColor(29); hym->Fit("gaus");
 259   c1->cd(3); hy3->Draw();
 260   c1->cd(4); hy2->Draw();
 261
 262 }
 263
 264
 265
 266
 267
 268