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 }