]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/Macros/AliTRDclusterErrors.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Move AliTRDclusterError.h into macro
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDclusterErrors.C
CommitLineData
d8d3b5b8 1#if !defined(__CINT__) || defined(__MAKECINT__)
aec95541 2#include <Riostream.h>
59dfcadd 3
d8d3b5b8 4#include "AliTRDtracker.h"
5#include "AliTRDclusterMI.h"
6#include "AliTRDhit.h"
7#include "AliTRDv1.h"
8#include "AliTRDgeometry.h"
9#include "AliTRDgeometryDetail.h"
10#include "AliTRDparameter.h"
aec95541 11#include "AliTRDclusterCorrection.h"
d8d3b5b8 12#include "alles.h"
13#include "TFile.h"
14#include "TStopwatch.h"
15#include "Rtypes.h"
16#include "TTree.h"
aec95541 17
289478c7 18#include "AliRunLoader.h"
19#include "AliStack.h"
289478c7 20#include "TF1.h"
21#include "AliTrackReference.h"
d8d3b5b8 22#endif
23#include "AliTRDclusterErrors.h"
289478c7 24
25
88cadd17 26AliTRDclusterCorrection * gCorrection;
27
28void ReadCorrection(){
29 TFile f("TRDcorrection.root");
30 gCorrection= (AliTRDclusterCorrection *)f.Get("TRDcorrection");
31 if (gCorrection==0){
32 printf("Correction not found");
33 }
34}
35
36
37class AliTRDExactPoint: public TObject {
38 public :
39 AliTRDExactPoint();
40 Float_t fTX; //x in rotated coordinate in the center of time bin
41 Float_t fTY; //y
42 Float_t fTZ; //z
43 Float_t fTAY; //angle y
44 Float_t fTAZ; //angle z
45 Float_t fGx;
46 Float_t fGy;
47 Float_t fGz;
48 //
49 void SetReference(AliTrackReference *ref);
50 Float_t fTRefAngleY;
51 Float_t fRefPos[3];
52 Float_t fRefMom[3];
53 //
54 Int_t fDetector; // detector (chamber)
55 Int_t fLocalTimeBin; // local time bin
56 Int_t fPlane; // plane (layer)
57 Int_t fSector; // segment
58 Int_t fPlaneMI;
59 //
60 Float_t fTQ;
61 Float_t fTPrim;
62 //
63 ClassDef(AliTRDExactPoint,1)
64};
65
66class AliTRDCI: public TObject {
67 public :
68 AliTRDCI(){;}
69 virtual ~AliTRDCI(){;}
70 //
71 AliTRDclusterMI fCl;
72 AliTRDExactPoint fEp;
73 TParticle fP;
74 Char_t fStatus;
75 Int_t fNClusters;
76 //
77 Float_t fDYtilt;
78 Float_t fTDistZ;
79 //
80 Int_t fNTracks;
81 Float_t fPt;
82 Float_t fCharge;
83 Bool_t fIsPrim;
84 Float_t fCorrection;
85 void Update();
86 ClassDef(AliTRDCI,1)
87};
88
89class AliTRDClusterErrAnal: public TObject{
90public:
91 AliTRDClusterErrAnal(Char_t *chloader = "galice.root");
92 void SetIO(Int_t event);
93 Int_t Analyze(Int_t trackmax);
94 void LoadClusters();
95 void MakeExactPoints(Int_t trackmax);
96 void SortReferences();
97 AliTrackReference * FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax=10.);
98public:
99 AliRunLoader * fRunLoader;
100 AliLoader * fTRDLoader;
101 AliTRDparameter *fParam;
102 AliTRDgeometry *fGeometry;
103 TTree * fHitTree;
104 TTree * fClusterTree;
105 TTree * fReferenceTree;
106 AliTRDv1 * fTRD;
107 //
108 TTree * fTreeA;
109 TFile * fFileA;
110 AliTRDtracker *fTracker;
111 AliStack *fStack;
112 TObjArray fClusters[12][100][18]; //first plane, second time bin
113 TObjArray fExactPoints;
114 TObjArray *fReferences;
115
116 ClassDef(AliTRDClusterErrAnal,1)
117};
118
119
120class AliTRDClusterErrDraw: public TObject{
121public:
122 TTree * fTree;
123 AliTRDclusterCorrection* MakeCorrection(TTree * tree, Float_t offset);
124 static TH1F * ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin=10, Float_t ampmax=300);
125 static TH1F * ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min=-0.5, Float_t max=0.5);
126 static TH1F * ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min=-1., Float_t max=1.);
127 static void AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle);
128 static TH1F* CreateEffHisto(TH1F* hGen, TH1F* hRec);
129 static TH1F* CreateResHisto(TH2F* hRes2, Bool_t draw = kTRUE, Bool_t drawBinFits = kTRUE,
130 Bool_t overflowBinFits = kFALSE);
131 ClassDef(AliTRDClusterErrDraw,1)
132};
133
134
135
136
137AliTRDExactPoint::AliTRDExactPoint()
138{
139 fTX=fTY=fTZ=fTAZ=fTAY=fGx=fGy=fGz=fTRefAngleY=0;
140 fRefPos[0]=fRefPos[1]=fRefPos[2]=fRefMom[0]=fRefMom[1]=fRefMom[2]=0;
141 fDetector=fLocalTimeBin=fPlane=fSector=fPlaneMI=0;
142 fTQ=fTPrim=0;
143}
144
145void AliTRDExactPoint::SetReference(AliTrackReference *ref){
146 fRefPos[0] = ref->X();
147 fRefPos[1] = ref->Y();
148 fRefPos[2] = ref->Z();
149 //
150 fRefMom[0] = ref->Px();
151 fRefMom[1] = ref->Py();
152 fRefMom[2] = ref->Pz();
153}
154
155
156void AliTRDCI::Update()
157{
158 //
159 //thanks to root
160 fPt = fP.Pt();
161 fCharge = fP.GetPDG()->Charge();
162 fIsPrim = (fP.GetMother(0)>0)? kFALSE :kTRUE;
163 fCorrection = gCorrection->GetCorrection(fEp.fPlane,fCl.fTimeBin,fEp.fTAY);
164}
165
166
167/*
168//example seesion
169
170.L AliGenInfo.C+
171.L AliTRDclusterErrors.C+
172gCorrection = AliTRDclusterCorrection::GetCorrection();
173AliTRDClusterErrAnal ana;
174ana.Analyze(10000000)
175
176
177
178.L AliGenInfo.C+
179.L AliTRDclusterErrors.C+
180TFile f("trdclusteranal.root");
181TTree* tree = (TTree*)f.Get("trdcl");
182AliComparisonDraw comp;
183comp->fTree = tree;
184
185
186tree->SetAlias("shapef","(1.-(0.8+0.06*(6-fEp.fPlane))*(fCl.fSigmaY2/(0.17+0.027*abs(fEp.fTAY))))");
187tree->SetAlias("shapes","0.08+0.3/sqrt(fCl.fQ)");
188tree->SetAlias("sfactor","shapef/shapes");
189
190tree->SetAlias("shapen","(fCl.fNPads-2.7-0.9*abs(fEp.fTAY))");
191
192tree->SetAlias("gshape","sfactor>-2&&fCl.fNPads<6&&shapen<1");
193
194
195tree->SetAlias("dy" , "fEp.fTY-fCl.fY-fDYtilt");
196tree->SetAlias("angle","abs(fEp.fTAY)");
197TCut cbase("cbase","(abs(fP.fPdgCode)!=11||fPt>0.2)&&fPt>0.1&&angle<2");
198
199tree->SetAlias("erry0","(0.028+0.07*angle)");
200tree->SetAlias("erry1","erry0*(0.9+15./fCl.fQ)");
201tree->SetAlias("erry2","erry1*(0.8+0.5*max(-sfactor,0))");
202
203
204
205
206TH1F his("resy","resy",100,-0.2,0.2);
207comp->fTree->Draw("dy:0.028*fEp.fTAY","fStatus==0 && abs(dy)<1.0&&fNTracks<2&&angle<2&&gshape","")
208comp->DrawXY("sqrt(fCl.fQ)","dy","fStatus==0"+cbase,"gshape",10,0,20,-0.7,0.7);
209
210comp->DrawXY("angle","dy/erry1","fStatus==0"+cbase,"gshape",10,0,2,-5.7,5.7);
211comp->DrawXY("sqrt(cl->fQ)","dy/err1","fStatus==0"+cbase,"gshape",10,2,20,-5.7,5.7);
212
213
214
215AliTRDClusterErrDraw::ResDyVsAmp(tree,"abs(dy)<0.4&&abs(fEp.fTAY)<0.2&&fNTracks<2&&fPt>0.1"+cbase,-0.03);
216AliTRDClusterErrDraw::ResDyVsRelPos(tree,"abs(dy)<0.4&&abs(fEp.fTAY)<0.2&&fNTracks<2&&fPt>0.2"+cbase,-0.03);
217AliTRDClusterErrDraw::ResDyVsAngleY(tree,"abs(fEp.fTY-fCl.fY+fDYtilt)<0.4",0.);
218
219AliTRDclusterCorrection * cor = AliTRDClusterErrDraw::MakeCorrection(tree,0.134);
220tree->Draw("sqrt(fCl.fRmsY)","fStatus==0&&abs(fEp.fTY-fCl.fY+fDYtilt)<0.4&&abs(fEp.fTAY)<0.2&&fNTracks<2&&fPt>0.2&&fCl.fNPads==0");
221
222tree->Draw("fEp.fTY-fCl.fY+fDYtilt:fCl.fTimeBin","fStatus==0&&abs(fEp.fTY-fCl.fY+fDYtilt)<0.5&&abs(fEp.fTAY)<0.3&&fEp.fTAY>0.13","prof")
223
224
225*/
226
227
289478c7 228ClassImp(AliTRDCI)
229ClassImp(AliTRDExactPoint)
230ClassImp(AliTRDClusterErrAnal)
231ClassImp(AliTRDClusterErrDraw)
232
233
234
235AliTRDClusterErrAnal::AliTRDClusterErrAnal(Char_t *chloader )
236{
237 //
238 //SET Input loaders
239 if (gAlice){
240 delete gAlice->GetRunLoader();
241 delete gAlice;
242 gAlice = 0x0;
243 }
244 fRunLoader = AliRunLoader::Open(chloader);
245 if (fRunLoader == 0x0){
246 cerr<<"Can not open session"<<endl;
247 return;
248 }
249 fTRDLoader = fRunLoader->GetLoader("TRDLoader");
250 if (fTRDLoader == 0x0){
251 cerr<<"Can not get TRD Loader"<<endl;
252 return ;
253 }
254 if (fRunLoader->LoadgAlice()){
255 cerr<<"Error occured while l"<<endl;
256 return;
257 }
258 fRunLoader->CdGAFile();
259 fTracker = new AliTRDtracker(gFile);
260 fParam = (AliTRDparameter*) gFile->Get("TRDparameter");
261 fGeometry = new AliTRDgeometryDetail();
262 fTRD = (AliTRDv1*) gAlice->GetDetector("TRD");
263 //
264 AliTRDCI * clinfo = new AliTRDCI();
265 fFileA = new TFile("trdclusteranal.root","recreate");
266 fFileA->cd();
267 fTreeA = new TTree("trdcl","trdcl");
268 fTreeA->Branch("trdcl","AliTRDCI",&clinfo);
269 delete clinfo;
270}
59dfcadd 271
289478c7 272void AliTRDClusterErrAnal::SetIO(Int_t event)
273{
274 //
275 //set input output for given event
276 fRunLoader->SetEventNumber(event);
277 fRunLoader->LoadHeader();
278 fRunLoader->LoadKinematics();
279 fRunLoader->LoadTrackRefs();
280 fTRDLoader->LoadHits();
281 fTRDLoader->LoadRecPoints("read");
282 //
283 fStack = fRunLoader->Stack();
284 fHitTree = fTRDLoader->TreeH();
285 fClusterTree = fTRDLoader->TreeR();
286 fReferenceTree = fRunLoader->TreeTR();
f007cb5f 287 fTracker->LoadClusters(fClusterTree);
289478c7 288 //
289}
59dfcadd 290
289478c7 291void AliTRDClusterErrAnal::LoadClusters()
292{
293 //
294 //
295 // Load clusters
296 TObjArray *ClusterArray = new TObjArray(400);
59dfcadd 297 TObjArray carray(2000);
289478c7 298 TBranch *branch=fClusterTree->GetBranch("TRDcluster");
299 if (!branch) {
300 Error("ReadClusters","Can't get the branch !");
301 return;
302 }
f007cb5f 303 Int_t over5 =0;
304 Int_t over10=0;
305
289478c7 306 branch->SetAddress(&ClusterArray);
307 Int_t nentries = (Int_t)fClusterTree->GetEntries();
308 for (Int_t i=0;i<nentries;i++){
309 fClusterTree->GetEvent(i);
310 Int_t nCluster = ClusterArray->GetEntriesFast();
311 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
312 AliTRDcluster* c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
313 carray.AddLast(c);
314 ClusterArray->RemoveAt(iCluster);
f007cb5f 315 if (c->GetQ()>5) over5++;
316 if (c->GetQ()>10) over10++;
289478c7 317 }
318 }
59dfcadd 319 Int_t nClusters = carray.GetEntriesFast();
f007cb5f 320 printf("Total number of clusters %d\t%d\t%d\n", nClusters,over5,over10);
289478c7 321 //
322 //
323 //SORT clusters
324 //
325 Int_t all=0;
326 for (Int_t i=0;i<nClusters;i++){
327 AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);
328 Int_t plane = fGeometry->GetPlane(cl->GetDetector());
329 if (plane>=12) continue;
330 Int_t time = cl->GetLocalTimeBin();
331 if (time>=100) continue;
332 Int_t sector = fGeometry->GetSector(cl->GetDetector());
333 if (sector>=18){
334 printf("problem1\n");
335 continue;
336 }
337 fClusters[plane][time][sector].AddLast(cl);
338 all++;
59dfcadd 339 }
289478c7 340}
59dfcadd 341
289478c7 342void AliTRDClusterErrAnal::SortReferences()
343{
344 //
345 //
346 //
347 printf("Sorting references\n");
348 fReferences = new TObjArray;
349 Int_t ntracks = fStack->GetNtrack();
350 fReferences->Expand(ntracks);
351 Int_t nentries = (Int_t)fReferenceTree->GetEntries();
352 TClonesArray * arr = new TClonesArray("AliTrackReference");
353 TBranch * br = fReferenceTree->GetBranch("TRD");
354 br->SetAddress(&arr);
355 //
356 TClonesArray *labarr=0;
aec95541 357 Int_t nreferences=0;
358 Int_t nreftracks=0;
289478c7 359 for (Int_t iprim=0;iprim<nentries;iprim++){
360 if (br->GetEntry(iprim)){
361 for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
362 AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
363 if (!ref) continue;
364 Int_t lab = ref->GetTrack();
365 if ( (lab<0) || (lab>ntracks)) continue;
366 //
367 if (fReferences->At(lab)==0) {
368 labarr = new TClonesArray("AliTrackReference");
369 fReferences->AddAt(labarr,lab);
aec95541 370 nreftracks++;
289478c7 371 }
372 TClonesArray &larr = *labarr;
373 new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
aec95541 374 nreferences++;
289478c7 375 }
376 }
377 }
aec95541 378 printf("Total number of references = \t%d\n", nreferences);
379 printf("Total number of tracks with references = \t%d\n", nreftracks);
380 printf("End - Sorting references\n");
381
289478c7 382}
59dfcadd 383
289478c7 384AliTrackReference * AliTRDClusterErrAnal::FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax)
385{
386 //
387 //
388 //
389 if (fReferences->At(lab)==0) return 0;
390 AliTrackReference *nearest=0;
391 TClonesArray * arr = (TClonesArray *)fReferences->At(lab);
392 for (Int_t iref =0;iref<arr->GetEntriesFast();iref++){
393 AliTrackReference * ref = ( AliTrackReference *)arr->UncheckedAt(iref);
394 if (!ref) continue;
395 Float_t delta = (pos[0]-ref->X())*(pos[0]-ref->X());
396 delta += (pos[1]-ref->Y())*(pos[1]-ref->Y());
397 delta += (pos[2]-ref->Z())*(pos[2]-ref->Z());
398 delta = TMath::Sqrt(delta);
399 if (delta<dmax){
400 dmax=delta;
401 nearest = ref;
402 }
403 }
404 return nearest;
405}
59dfcadd 406
289478c7 407void AliTRDClusterErrAnal::MakeExactPoints(Int_t trackmax)
408{
409 //
410 //make exact points:-)
411 //
412 //
413
414 fExactPoints.Delete();
415 fExactPoints.Expand(fStack->GetNtrack());
416 //
417 Double_t fSum=0;
418 Double_t fSumQ =0;
419 Double_t fSumX=0;
420 Double_t fSumX2=0;
421 Double_t fSumXY=0;
422 Double_t fSumXZ=0;
423 Double_t fSumY=0;
424 Double_t fSumZ=0;
425 //
426 Int_t entries = Int_t(fHitTree->GetEntries());
427 printf("Number of primary entries\t%d\n",entries);
428 entries = TMath::Min(trackmax,entries);
429 Int_t nallpoints = 0;
430
431 Int_t nalltracks =0;
432 Int_t pointspertrack =0;
433
434 for (Int_t entry=0;entry<entries; entry++){
435 gAlice->ResetHits();
436 fHitTree->GetEvent(entry);
437 Int_t lastlabel = -1;
438 Int_t lastdetector = -1;
439 Int_t lasttimebin = -1;
440 Float_t lastpos[3];
441 //
442 for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); hit;
443 hit = (AliTRDhit *) fTRD->NextHit()) {
444 //
445 Int_t label = hit->Track();
446 TParticle * particle = fStack->Particle(label);
447 if (!particle) continue;
448 if (particle->Pt()<0.05) continue;
449 Int_t detector = hit->GetDetector();
450 Int_t plane = fGeometry->GetPlane(detector);
451 //
452 //
453 if (hit->GetCharge()==0) continue;
454 Float_t pos[3] = {hit->X(),hit->Y(),hit->Z()};
455 Int_t indexes[3];
456 fGeometry->Global2Detector(pos,indexes,fParam);
457 //
458 Float_t rot[3];
459 fGeometry->Rotate(detector,pos,rot);
aec95541 460 //rot[0] *=-1;
461 // rot[1] *=-1;
289478c7 462 //
463 //
464 Float_t time0 = fParam->GetTime0(plane);
465 Int_t timebin = Int_t(TMath::Nint(((time0 - rot[0])/fParam->GetTimeBinSize())+ fParam->GetTimeBefore())+0.1);
466 if (timebin<0) continue;
467 //
468 //
469 if (label!=lastlabel || detector != lastdetector || lasttimebin !=timebin){
470 //
471 if (label!=lastlabel){
472 fExactPoints.AddAt(new TClonesArray("AliTRDExactPoint",0),label);
473 //printf("new particle\t%d\n",label);
474 nalltracks++;
475 // printf("particle\t%d- hits\t%d\n",lastlabel, pointspertrack);
476 pointspertrack=0;
477 }
478
479 if ( (fSum>1) && lasttimebin>=0 && lasttimebin<fParam->GetTimeMax() ){
480 //if we have enough info for given layer time bin - store it
481 AliTrackReference * ref = FindNearestReference(lastlabel,lastpos,4.);
482 Float_t rotmom[3];
483 Float_t rotpos[3];
484 Float_t refangle=0;
485 if (ref){
486 Float_t mom[3] = {ref->Px(),ref->Py(),ref->Pz()};
487 Float_t refpos[3] = {ref->X(),ref->Y(),ref->Z()};
488 fGeometry->Rotate(detector,mom,rotmom);
489 fGeometry->Rotate(detector,refpos,rotpos);
490 refangle = rotmom[1]/rotmom[0];
491
492 }
493
494 Double_t ay,by,az,bz;
495 Double_t det = fSum*fSumX2-fSumX*fSumX;
496 if (TMath::Abs(det)> 0.000000000000001) {
497 by = (fSum*fSumXY-fSumX*fSumY)/det;
498 ay = (fSumX2*fSumY-fSumX*fSumXY)/det;
499
500 }else{
501 ay =fSumXY/fSumX;
502 by =0;
503 }
504 if (TMath::Abs(det)> 0.000000000000001) {
505 bz = (fSum*fSumXZ-fSumX*fSumZ)/det;
506 az = (fSumX2*fSumZ-fSumX*fSumXZ)/det;
507 }else{
508 az =fSumXZ/fSumX;
509 bz =0;
510 }
511 //
512 Float_t lastplane = fGeometry->GetPlane(lastdetector);
513 Float_t time0 = fParam->GetTime0(lastplane);
aec95541 514 Float_t xcenter0 = time0 - (lasttimebin - fParam->GetTimeBefore()+0.5)*fParam->GetTimeBinSize();
515 Float_t xcenter = fTracker->GetX(0,lastplane,lasttimebin);
516 if (TMath::Abs(xcenter-xcenter0)>0.001){
517 printf("problem");
518 }
519
289478c7 520 Float_t ty = ay + by * xcenter;
521 Float_t tz = az + bz * xcenter;
522 //
523
524 TClonesArray * arr = (TClonesArray *) fExactPoints.At(label);
525 TClonesArray & larr= *arr;
526 Int_t arrent = arr->GetEntriesFast();
527 AliTRDExactPoint * point = new (larr[arrent]) AliTRDExactPoint;
528 nallpoints++;
529
530 if (ref){
531 point->SetReference(ref);
532 point->fTRefAngleY = rotmom[1]/rotmom[0];
533 }
534 point->fTX = xcenter;
535 point->fTY = ty;
536 point->fTZ = tz;
537 point->fTAY = by;
538 point->fTAZ = bz;
539 //
540 point->fGx = lastpos[0];
541 point->fGy = lastpos[1];
542 point->fGz = lastpos[2];
543
544 //
545 point->fDetector = lastdetector;
546 point->fLocalTimeBin = lasttimebin;
547 point->fPlane = fGeometry->GetPlane(lastdetector);
548 point->fSector = fGeometry->GetSector(lastdetector);
549 point->fPlaneMI = indexes[0];
550 //
551 point->fTPrim = fSum;
552 point->fTQ = fSumQ;
553 //
554 }
555 lastdetector = detector;
556 lastlabel = label;
557 lasttimebin = timebin;
558 fSum=fSumQ=fSumX=fSumX2=fSumXY=fSumXZ=fSumY=fSumZ=0.;
559 }
560 //
561 lastpos[0] = hit->X();
562 lastpos[1] = hit->Y();
563 lastpos[2] = hit->Z();
564 fSum++;
aec95541 565 fSumQ +=hit->GetCharge();
289478c7 566 fSumX +=rot[0];
567 fSumX2 +=rot[0]*rot[0];
568 fSumXY +=rot[0]*rot[1];
569 fSumXZ +=rot[0]*rot[2];
570 fSumY +=rot[1];
571 fSumZ +=rot[2];
572 pointspertrack++;
573 }
574 }
575 //
576 printf("Found %d exact points\n",nallpoints);
577}
59dfcadd 578
59dfcadd 579
59dfcadd 580
59dfcadd 581
59dfcadd 582
59dfcadd 583
289478c7 584Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {
59dfcadd 585
289478c7 586 //
587 // comparison works with both cluster types MI and old also
588 //dummy cluster to be fill if not cluster info
589 AliTRDclusterMI clmi;
590 TClass * classmi = clmi.IsA();
591 //
592 //SetOutput
593 AliTRDCI * clinfo = new AliTRDCI();
594 TBranch * clbr = fTreeA->GetBranch("trdcl");
595 clbr->SetAddress(&clinfo);
59dfcadd 596
289478c7 597 SetIO(0);
598 SortReferences();
599 MakeExactPoints(trackmax);
600 LoadClusters();
601 //
602 trackmax = fStack->GetNtrack();
603 //
604 // Get the number of entries in the hit tree
605 // (Number of primary particles creating a hit somewhere)
606 Int_t nTrack = (Int_t)fExactPoints.GetEntries();
aec95541 607 printf("Found %d charged in TRD in first %d particles", nTrack, trackmax);
289478c7 608 //
609
610 for (Int_t itrack = 0; itrack<trackmax; itrack++){
611 TClonesArray *arrpoints = (TClonesArray*)fExactPoints.At(itrack);
612
613 if (!arrpoints) continue;
614 //printf("new particle\t%d\n",itrack);
615 TParticle * particle = fStack->Particle(itrack);
616 if (!particle) continue;
617 //printf("founded in kine tree \t%d\n",itrack);
618 Int_t npoints = arrpoints->GetEntriesFast();
619 if (npoints<10) continue;
620 //printf("have enough points \t%d\t%d\n",itrack,npoints);
621
622 for (Int_t ipoint=0;ipoint<npoints;ipoint++){
623 AliTRDExactPoint * point = (AliTRDExactPoint*)arrpoints->UncheckedAt(ipoint);
624 if (!point) continue;
625 //
626 Int_t sec = fGeometry->GetSector(point->fDetector);
627 if (sec>18){
628 printf("problem2\n");
629 }
630 TObjArray & cllocal = fClusters[point->fPlane][point->fLocalTimeBin][sec];
631 Int_t nclusters = cllocal.GetEntriesFast();
632 Float_t maxdist = 10;
633 AliTRDcluster * nearestcluster =0;
f007cb5f 634 clinfo->fNClusters=0;
289478c7 635 //find nearest cluster to hit with given label
636 for (Int_t icluster =0; icluster<nclusters; icluster++){
637 AliTRDcluster * cluster = (AliTRDcluster*)cllocal.UncheckedAt(icluster);
638 if (!cluster) continue;
639 if ( (cluster->GetLabel(0)!= itrack) && (cluster->GetLabel(1)!= itrack)&&(cluster->GetLabel(2)!= itrack))
640 continue;
641 Float_t dist = TMath::Abs(cluster->GetY()-point->fTY);
f007cb5f 642 if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5 || dist>3.) continue;
643 clinfo->fNClusters++;
289478c7 644 if (dist<maxdist){
645 maxdist = dist;
646 nearestcluster = cluster;
647 }
648 }
649 //
650 clinfo->fEp = *point;
651 clinfo->fP = *particle;
652 if (!nearestcluster) {
653 clinfo->fStatus=1;
654 clinfo->fCl = clmi;
655 }
656 else{
657 clinfo->fStatus=0;
658 if (nearestcluster->IsA()==classmi){
659 clinfo->fCl =*((AliTRDclusterMI*)nearestcluster);
660 }
661 else{
662 clinfo->fCl = *nearestcluster;
59dfcadd 663 }
289478c7 664 //
665 Float_t dz = clinfo->fCl.GetZ()-point->fTZ;
666 Double_t h01 = sin(TMath::Pi() / 180.0 * fParam->GetTiltingAngle());
667 clinfo->fTDistZ = dz;
668 clinfo->fDYtilt = h01*dz*((point->fPlane%2)*2.-1.);
669 //
670 clinfo->fNTracks =1;
671 if (nearestcluster->GetLabel(1)>=0) clinfo->fNTracks++;
672 if (nearestcluster->GetLabel(2)>=0) clinfo->fNTracks++;
673 clinfo->Update();
674 }
675 //
676 fTreeA->Fill();
677 }
678 }
679
680
681 fFileA->cd();
682 fTreeA->Write();
683 fFileA->Close();
684 return 0;
685}
59dfcadd 686
aec95541 687AliTRDclusterCorrection* AliTRDClusterErrDraw::MakeCorrection(TTree * tree, Float_t offset)
688{
689 //
690 //
691 // make corrections
692 AliTRDclusterCorrection * cor = new AliTRDclusterCorrection;
693 cor->SetOffsetAngle(offset);
694 for (Int_t iplane=0;iplane<6;iplane++)
695 for (Int_t itime=0;itime<15;itime++)
696 for (Int_t iangle=0; iangle<20;iangle++){
697 Float_t angle = cor->GetAngle(iangle);
698 TH1F delta("delta","delta",30,-0.3,0.3);
699 char selection[100]="fStatus==0&&fNTracks<2";
700 char selectionall[1000];
701 sprintf(selectionall,"%s&&abs(fEp.fTAY-%f)<0.2&&fEp.fPlane==%d&&fCl.fTimeBin==%d&&fCl.fQ>20",
702 selection,angle,iplane,itime);
703 printf("\n%s",selectionall);
704 tree->Draw("fEp.fTY-fCl.fY+fDYtilt>>delta",selectionall);
705 gPad->Update();
706 printf("\nplane\t%d\ttime%d\tangle%f",iplane,itime,angle);
707 printf("\tentries%f\tmean\t%f\t%f",delta.GetEntries(),delta.GetMean(),delta.GetRMS());
708 cor->SetCorrection(iplane,itime,angle,delta.GetMean(),delta.GetRMS());
709 }
710 TFile * f = new TFile("TRDcorrection.root","new");
711 if (!f) f = new TFile("TRDcorrection.root","recreate");
712 f->cd();
713 cor->Write("TRDcorrection");
714 f->Close();
715 return cor;
716}
59dfcadd 717
289478c7 718TH1F * AliTRDClusterErrDraw::ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin, Float_t ampmax)
719{
720 //
721 //
722 TH2F hisdy("resy","resy",10,ampmin,ampmax,30,-0.3,0.3);
723 char expression[1000];
724 sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fQ>>resy",t0);
725 char selectionall[1000];
726 sprintf(selectionall,"fStatus==0&&%s",selection);
727 printf("%s\n",expression);
728 printf("%s\n",selectionall);
729 tree->Draw(expression,selectionall);
730 return CreateResHisto(&hisdy);
731}
59dfcadd 732
59dfcadd 733
289478c7 734TH1F * AliTRDClusterErrDraw::ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
735{
736 //
737 //
738 min *=128;
739 max *=128;
740 TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
741 char expression[1000];
742 sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fRelPos>>resy",t0);
743 char selectionall[1000];
744 sprintf(selectionall,"fStatus==0&&%s",selection);
745 printf("%s\n",expression);
746 printf("%s\n",selectionall);
747 tree->Draw(expression,selectionall);
748 return CreateResHisto(&hisdy);
749}
59dfcadd 750
751
289478c7 752TH1F * AliTRDClusterErrDraw::ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
753{
754 //
755 //
756 TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
59dfcadd 757
289478c7 758 char expression[1000];
759 sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%f*fEp.fTAY:fEp.fTAY>>resy",t0);
760 char selectionall[1000];
761 sprintf(selectionall,"fStatus==0&&%s",selection);
59dfcadd 762
289478c7 763 tree->Draw(expression,selectionall);
764 return CreateResHisto(&hisdy);
765}
59dfcadd 766
289478c7 767void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
768{
769 histo->GetXaxis()->SetTitle(xAxisTitle);
770 histo->GetYaxis()->SetTitle(yAxisTitle);
771}
59dfcadd 772
59dfcadd 773
289478c7 774TH1F* AliTRDClusterErrDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
775{
776 Int_t nBins = hGen->GetNbinsX();
777 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
778 hEff->SetTitle("");
779 hEff->SetStats(kFALSE);
780 hEff->SetMinimum(0.);
781 hEff->SetMaximum(110.);
782
783 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
784 Double_t nGen = hGen->GetBinContent(iBin);
785 Double_t nRec = hRec->GetBinContent(iBin);
786 if (nGen > 0) {
787 Double_t eff = nRec/nGen;
788 hEff->SetBinContent(iBin, 100. * eff);
789 Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
790 // if (error == 0) error = sqrt(0.1/nGen);
791 //
792 if (error == 0) error = 0.0001;
793 hEff->SetBinError(iBin, 100. * error);
794 } else {
795 hEff->SetBinContent(iBin, 100. * 0.5);
796 hEff->SetBinError(iBin, 100. * 0.5);
797 }
798 }
59dfcadd 799
289478c7 800 return hEff;
801}
59dfcadd 802
803
59dfcadd 804
289478c7 805TH1F* AliTRDClusterErrDraw::CreateResHisto(TH2F* hRes2, Bool_t draw, Bool_t drawBinFits,
806 Bool_t overflowBinFits)
807{
808 TVirtualPad* currentPad = gPad;
809 TAxis* axis = hRes2->GetXaxis();
810 Int_t nBins = axis->GetNbins();
811 TH1F* hRes, *hMean;
812 if (axis->GetXbins()->GetSize()){
813 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
814 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
815 }
816 else{
817 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
818 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
59dfcadd 819
289478c7 820 }
821 hRes->SetStats(false);
822 hRes->SetOption("E");
823 hRes->SetMinimum(0.);
824 //
825 hMean->SetStats(false);
826 hMean->SetOption("E");
827
828 // create the fit function
829 //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
830 TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
831
832 fitFunc->SetLineWidth(2);
833 fitFunc->SetFillStyle(0);
834 // create canvas for fits
835 TCanvas* canBinFits = NULL;
836 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
837 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
838 Int_t ny = (nPads-1) / nx + 1;
839 if (drawBinFits) {
840 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
841 if (canBinFits) delete canBinFits;
842 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
843 canBinFits->Divide(nx, ny);
844 }
59dfcadd 845
289478c7 846 // loop over x bins and fit projection
847 Int_t dBin = ((overflowBinFits) ? 1 : 0);
848 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
849 if (drawBinFits) canBinFits->cd(bin + dBin);
850 TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
851 //
852 if (hBin->GetEntries() > 10) {
853 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS(),0.02*hBin->GetMaximum());
854 hBin->Fit(fitFunc,"s");
855 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
856
857 if (sigma > 0.){
858 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
859 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
59dfcadd 860 }
289478c7 861 else{
862 hRes->SetBinContent(bin, 0.);
863 hMean->SetBinContent(bin,0);
864 }
865 hRes->SetBinError(bin, fitFunc->GetParError(2));
866 hMean->SetBinError(bin, fitFunc->GetParError(1));
867
868 //
869 //
870
871 } else {
872 hRes->SetBinContent(bin, 0.);
873 hRes->SetBinError(bin, 0.);
874 hMean->SetBinContent(bin, 0.);
875 hMean->SetBinError(bin, 0.);
59dfcadd 876 }
289478c7 877
878
879 if (drawBinFits) {
880 char name[256];
881 if (bin == 0) {
882 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
883 } else if (bin == nBins+1) {
884 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
885 } else {
886 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
887 axis->GetTitle(), axis->GetBinUpEdge(bin));
888 }
889 canBinFits->cd(bin + dBin);
890 hBin->SetTitle(name);
891 hBin->SetStats(kTRUE);
892 hBin->DrawCopy("E");
893 canBinFits->Update();
894 canBinFits->Modified();
895 canBinFits->Update();
896 }
897
898 delete hBin;
899 }
900
901 delete fitFunc;
902 currentPad->cd();
903 if (draw){
904 currentPad->Divide(1,2);
905 currentPad->cd(1);
906 hRes->Draw();
907 currentPad->cd(2);
908 hMean->Draw();
59dfcadd 909 }
910
289478c7 911 return hRes;
59dfcadd 912}