4 #include "AliTRDtracker.h"
5 #include "AliTRDcluster.h"
8 #include "AliTRDgeometry.h"
9 #include "AliTRDparameter.h"
13 #include "TStopwatch.h"
17 void AliTRDclusterErrors() {
19 Int_t No_of_tracks_to_analyze = 40;
21 const Int_t tbpp = 15;
22 const Int_t nPlanes = 6;
23 const Int_t ntb = tbpp * nPlanes;
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.);
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);
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");
47 // Dynamically link some shared libs
48 if (gClassTable->GetID("AliRun") < 0) {
49 gROOT->LoadMacro("loadlibs.C");
51 cout << "Loaded shared libraries" << endl;
56 Char_t *clusterfile = "AliTRDclusters.root";
59 TObjArray carray(2000);
60 TObjArray *ClustersArray = &carray;
61 TFile *geofile =TFile::Open("AliTRDclusters.root");
63 AliTRDparameter *par = (AliTRDparameter*) geofile->Get("TRDparameter");
64 AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry");
66 AliTRDtracker *Tracker = new AliTRDtracker(geofile);
67 Tracker->SetEventNumber(nEvent);
68 Tracker->ReadClusters(ClustersArray,clusterfile);
69 Int_t nClusters = carray.GetEntriesFast();
71 printf("Total number of clusters %d \n", nClusters);
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);
77 cout << "Open the ALIROOT-file " << alifile << endl;
78 gafl = new TFile(alifile);
81 cout << alifile << " is already open" << endl;
84 // Get AliRun object from file or create it if not on file
85 gAlice = (AliRun*) gafl->Get("gAlice");
87 cout << "AliRun object found on file" << endl;
89 gAlice = new AliRun("gAlice","Alice test program");
91 AliTRDv1 *fTRD = (AliTRDv1*) gAlice->GetDetector("TRD");
93 // Import the Trees for the event nEvent in the file
94 Int_t nparticles = gAlice->GetEvent(nEvent);
96 TObjArray *particles=gAlice->Particles();
98 TTree *hitTree = gAlice->TreeH();
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);
108 // Loop through particles and fill histoes
113 Float_t clusterY[ntb];
114 Float_t clusterZ[ntb];
115 Float_t clusterQ[ntb];
116 Float_t clusterSigmaY[ntb];
121 Int_t track_index[3];
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);
128 for(Int_t plane = 0; plane < nPlanes; plane++) {
129 for(Int_t ltb = 14; ltb > -1; ltb--) {
131 if(ii >= nTrack) continue;
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;
139 Int_t gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
140 Double_t x = Tracker->GetX(0,plane,ltb);
142 // loop through clusters
144 Bool_t cluster_found = kFALSE;
147 for (Int_t i = 0; i < nClusters; i++) {
149 AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);
151 nw = cl->GetDetector();
153 nw = fGeo->GetPlane(nw);
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;
161 nw = cl->GetLocalTimeBin();
162 if(nw != ltb) continue;
164 clusterY[gtb] = cl->GetY();
165 clusterZ[gtb] = cl->GetZ();
166 clusterQ[gtb] = cl->GetQ();
168 clusterSigmaY[gtb] = TMath::Sqrt(cl->GetSigmaY2());
169 cluster_found = kTRUE;
173 if(!cluster_found) continue;
177 // nBytes += hitTree->GetEvent(nPrimaries - ii - 1);
178 nBytes += hitTree->GetEvent(nTrack - ii - 1);
180 // Loop through the TRD hits
181 Bool_t found_hit = kFALSE;
184 for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
186 hit = (AliTRDhit *) fTRD->NextHit()) {
188 if(nw != ii) continue;
189 det = hit->GetDetector();
190 nw = fGeo->GetPlane(det);
191 if(nw != plane) continue;
197 fGeo->Rotate(det,pos,rot);
199 if(TMath::Abs(rot[0]-x) > 0.01) continue;
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;
211 if(!found_hit) continue;
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]);
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]);
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]);
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]);
243 gStyle->SetOptStat(1);
244 gStyle->SetOptFit(1);
246 TCanvas* c = new TCanvas("c", "c", 110, 110, 810, 840);
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();
254 TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940);
255 c1->SetFillColor(10);
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();