]> git.uio.no Git - u/mrichter/AliRoot.git/blob - TRD/Macros/AliTRDclusterErrors.C
A few mods in the GRP preprocessor (to get some greens finally)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDclusterErrors.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <Riostream.h>
3
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"
11 #include "AliTRDclusterCorrection.h"
12 #include "alles.h"
13 #include "TFile.h"
14 #include "TStopwatch.h"
15 #include "Rtypes.h"
16 #include "TTree.h"
17
18 #include "AliRunLoader.h"
19 #include "AliStack.h"
20 #include "TF1.h"
21 #include "AliTrackReference.h"
22 #endif    
23 #include "AliTRDclusterErrors.h"
24
25
26 AliTRDclusterCorrection * gCorrection;
27
28 void 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
37 class 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
66 class 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
89 class AliTRDClusterErrAnal: public TObject{
90 public: 
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.);
98 public:
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
120 class AliTRDClusterErrDraw: public TObject{
121 public:
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
137 AliTRDExactPoint::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
145 void 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
156 void 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+
172 gCorrection  = AliTRDclusterCorrection::GetCorrection();
173 AliTRDClusterErrAnal ana;
174 ana.Analyze(10000000)
175
176
177
178 .L AliGenInfo.C+
179 .L AliTRDclusterErrors.C+
180 TFile f("trdclusteranal.root");
181 TTree* tree = (TTree*)f.Get("trdcl");
182 AliComparisonDraw comp;
183 comp->fTree = tree;
184
185
186 tree->SetAlias("shapef","(1.-(0.8+0.06*(6-fEp.fPlane))*(fCl.fSigmaY2/(0.17+0.027*abs(fEp.fTAY))))");
187 tree->SetAlias("shapes","0.08+0.3/sqrt(fCl.fQ)");
188 tree->SetAlias("sfactor","shapef/shapes");
189
190 tree->SetAlias("shapen","(fCl.fNPads-2.7-0.9*abs(fEp.fTAY))");
191
192 tree->SetAlias("gshape","sfactor>-2&&fCl.fNPads<6&&shapen<1");
193
194
195 tree->SetAlias("dy"    , "fEp.fTY-fCl.fY-fDYtilt");
196 tree->SetAlias("angle","abs(fEp.fTAY)");
197 TCut cbase("cbase","(abs(fP.fPdgCode)!=11||fPt>0.2)&&fPt>0.1&&angle<2");
198
199 tree->SetAlias("erry0","(0.028+0.07*angle)");
200 tree->SetAlias("erry1","erry0*(0.9+15./fCl.fQ)");
201 tree->SetAlias("erry2","erry1*(0.8+0.5*max(-sfactor,0))"); 
202
203
204
205
206 TH1F his("resy","resy",100,-0.2,0.2);
207 comp->fTree->Draw("dy:0.028*fEp.fTAY","fStatus==0 && abs(dy)<1.0&&fNTracks<2&&angle<2&&gshape","")
208 comp->DrawXY("sqrt(fCl.fQ)","dy","fStatus==0"+cbase,"gshape",10,0,20,-0.7,0.7);
209
210 comp->DrawXY("angle","dy/erry1","fStatus==0"+cbase,"gshape",10,0,2,-5.7,5.7);
211 comp->DrawXY("sqrt(cl->fQ)","dy/err1","fStatus==0"+cbase,"gshape",10,2,20,-5.7,5.7);
212
213
214
215 AliTRDClusterErrDraw::ResDyVsAmp(tree,"abs(dy)<0.4&&abs(fEp.fTAY)<0.2&&fNTracks<2&&fPt>0.1"+cbase,-0.03);
216 AliTRDClusterErrDraw::ResDyVsRelPos(tree,"abs(dy)<0.4&&abs(fEp.fTAY)<0.2&&fNTracks<2&&fPt>0.2"+cbase,-0.03);
217 AliTRDClusterErrDraw::ResDyVsAngleY(tree,"abs(fEp.fTY-fCl.fY+fDYtilt)<0.4",0.);
218
219 AliTRDclusterCorrection * cor = AliTRDClusterErrDraw::MakeCorrection(tree,0.134);
220 tree->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
222 tree->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
228 ClassImp(AliTRDCI)
229 ClassImp(AliTRDExactPoint)
230 ClassImp(AliTRDClusterErrAnal)
231 ClassImp(AliTRDClusterErrDraw)
232
233
234
235 AliTRDClusterErrAnal::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 }
271
272 void 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();
287   fTracker->LoadClusters(fClusterTree);
288   //
289 }
290
291 void AliTRDClusterErrAnal::LoadClusters()
292 {
293   //
294   //
295   // Load clusters  
296   TObjArray *ClusterArray = new TObjArray(400);   
297   TObjArray carray(2000);
298   TBranch *branch=fClusterTree->GetBranch("TRDcluster");
299   if (!branch) {
300     Error("ReadClusters","Can't get the branch !");
301     return;
302   }
303   Int_t over5 =0;
304   Int_t over10=0;
305
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);
315       if (c->GetQ()>5)  over5++;
316       if (c->GetQ()>10) over10++;
317     }
318   }
319   Int_t nClusters = carray.GetEntriesFast();
320   printf("Total number of clusters %d\t%d\t%d\n", nClusters,over5,over10);
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++;
339   }
340 }
341
342 void 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;
357   Int_t nreferences=0;
358   Int_t nreftracks=0;
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);
370           nreftracks++;
371         }
372         TClonesArray &larr = *labarr;
373         new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
374         nreferences++;
375       }
376     }
377   }
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   
382 }
383
384 AliTrackReference * 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 }
406
407 void 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);
460       //rot[0] *=-1;
461       //  rot[1] *=-1;
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);
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           
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++;
565       fSumQ  +=hit->GetCharge();
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 }
578
579
580
581
582
583
584 Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {  
585
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);
596
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();
607   printf("Found %d charged in TRD in first %d particles", nTrack, trackmax);
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;
634       clinfo->fNClusters=0;
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);
642         if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5 || dist>3.) continue; 
643         clinfo->fNClusters++;
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;
663         }
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 }
686
687 AliTRDclusterCorrection*   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 }
717
718 TH1F * 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 }
732
733
734 TH1F * 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 }
750
751
752 TH1F * 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);
757
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);
762
763   tree->Draw(expression,selectionall);
764   return CreateResHisto(&hisdy);
765 }
766
767 void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
768 {
769   histo->GetXaxis()->SetTitle(xAxisTitle);
770   histo->GetYaxis()->SetTitle(yAxisTitle);
771 }
772
773
774 TH1F* 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   }
799
800   return hEff;
801 }
802
803
804
805 TH1F* 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());
819
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   }
845
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));    
860       }
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.);
876     }
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();
909   }
910   
911   return hRes;
912 }