Remove AliTRDclusterMI, AliTRDclusterCorrection, AliTRDmcTrack
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Feb 2008 09:19:17 +0000 (09:19 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Feb 2008 09:19:17 +0000 (09:19 +0000)
TRD/Macros/AliTRDbackTrackAnalysis.C [deleted file]
TRD/Macros/AliTRDclusterErrors.C [deleted file]
TRD/Macros/AliTRDsaveMCtracks.C [deleted file]
TRD/Macros/AliTRDsaveTrackableSeeds.C [deleted file]
TRD/Macros/AliTRDseedsAnalysis.C [deleted file]

diff --git a/TRD/Macros/AliTRDbackTrackAnalysis.C b/TRD/Macros/AliTRDbackTrackAnalysis.C
deleted file mode 100644 (file)
index 6d2ddf4..0000000
+++ /dev/null
@@ -1,606 +0,0 @@
-#ifndef __CINT__
-  #include <iostream.h>
-  #include "AliTRDtracker.h"
-  #include "AliTRDcluster.h" 
-  #include "AliTRDv1.h"
-  #include "AliTRDgeometry.h"    
-  #include "AliTRDparameter.h"    
-  #include "alles.h"  
-  #include "AliTRDmcTrack.h"
-  #include "AliTRDtrack.h"
-  
-  #include "TFile.h"
-  #include "TParticle.h"
-  #include "TStopwatch.h"
-#endif
-
-void AliTRDbackTrackAnalysis() {
-
-  //  const Int_t nPrimaries = 84210/400;
-  const Int_t nPrimaries = 84210/16;
-
-
-  Float_t Pt_min = 0.;
-  Float_t Pt_max = 20.;
-
-  TH1F *hp=new TH1F("hp","PHI resolution",100,-20.,20.); hp->SetFillColor(4);
-  TH1F *hl=new TH1F("hl","LAMBDA resolution",100,-100,100);hl->SetFillColor(4);
-  TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
-
-  TH1F *hpta=new TH1F("hpta","norm. Pt resolution",50,-5.,5.);
-  hpta->SetFillColor(17);
-  hpta->SetXTitle("(%)");
-
-  TH1F *hla=new TH1F("hla","norm. Lambda resolution",50,-5.,5.);
-  hla->SetFillColor(17);
-  TH1F *hya=new TH1F("hya","norm. Y resolution",50,-5.,5.);
-  hya->SetFillColor(17);
-  TH1F *hza=new TH1F("hza","norm. Z resolution",50,-5.,5.);
-  hza->SetFillColor(17);
-
-  hpt->SetFillColor(2);             
-
-  TH1F *hy=new TH1F("hy","Y resolution",50,-0.5,0.5);hy->SetLineColor(4);
-  hy->SetLineWidth(2);
-  hy->SetXTitle("(cm)");
-
-  TH1F *hz=new TH1F("hz","Z resolution",80,-4,4);hz->SetLineColor(2);
-  hz->SetLineWidth(2);
-  hz->SetXTitle("(cm)");
-
-  TH1F *hx=new TH1F("hx","X(out)",150,250.,400.); hx->SetFillColor(4);
-
-  const Int_t nPtSteps = 30;
-  const Float_t maxPt = 3.;
-
-  TH1F *hgood=new TH1F("hgood","long TRD tracks, TPC seeds",nPtSteps,0,maxPt);
-  hgood->SetYTitle("Counts");
-  hgood->SetXTitle("Pt (GeV/c)");
-  hgood->SetLineColor(3);
-
-  TH1F *hGoodAndSeed = new TH1F("GoodAndSeed","TPC seed and good",nPtSteps,0,maxPt);
-  hGoodAndSeed->SetLineColor(2);
-
-  TH2F *hRTvsMC = new TH2F("RTvsMC","RTvsMC",100,-4.5,95.5,100,-4.5,95.5); 
-  hRTvsMC->SetMarkerColor(4);
-  hRTvsMC->SetMarkerSize(2);
-
-  TH2F *hXvsMCX = new TH2F("XvsMCX","XvsMCX",150,250,400,150,250,400); 
-  hXvsMCX->SetMarkerColor(4);
-  hXvsMCX->SetMarkerSize(2);
-
-  TH2F *h2seed = new TH2F("2dGood","TPC seeds",60,0.,60,50,0.,3.); 
-  h2seed->SetMarkerColor(4);
-  h2seed->SetMarkerSize(2);
-
-  TH2F *h2lost = new TH2F("2dSeedAndGood","SeedButNotGood",60,0.,60.,50,0.,3.);  
-  h2lost->SetMarkerColor(2);
-  h2lost->SetMarkerSize(2);
-
-
-  TH1F *hseed=new TH1F("seed","TPC seeds",nPtSteps,0,maxPt);
-  hseed->SetLineColor(4);
-
-
-  TH1F *hfound=new TH1F("hfound","Found tracks",nPtSteps,0,maxPt);  
-  hfound->SetYTitle("Counts");
-  hfound->SetXTitle("Pt (GeV/c)");
-
-  TH1F *heff=new TH1F("heff","Matching Efficiency",nPtSteps,0,maxPt); // efficiency, good tracks
-  heff->SetLineColor(4); heff->SetLineWidth(2);  
-  heff->SetXTitle("Pt, GeV/c"); 
-  heff->SetYTitle("Efficiency"); 
-
-  TH1F *hSeedEff = new TH1F("hSeedEff","TPC Efficiency",nPtSteps,0,maxPt); 
-  hSeedEff->SetLineColor(4); hSeedEff->SetLineWidth(2);  
-  hSeedEff->SetXTitle("Pt, GeV/c"); 
-  hSeedEff->SetYTitle("Efficiency"); 
-
-
-  TH1F *hol=new TH1F("hol","Overlap fraction",105,-2.5,102.5); 
-  hol->SetLineColor(4); hol->SetLineWidth(2);  
-  hol->SetXTitle("Fraction,(%)"); 
-  hol->SetYTitle("Counts"); 
-
-  TH1F *hend=new TH1F("end","missing tail",80,-10.5,69.5);
-
-  TH1F *hFraction=new TH1F("fraction","Fraction of found clusters",110,0,1.1);
-  TH1F *hCorrect=new TH1F("correct","Fraction of correct clusters",110,0,1.1);
-
-
-  Int_t nEvent = 0;
-  const Int_t maxIndex = nPrimaries;
-  Bool_t seedLabel[maxIndex];
-  Int_t mcIndex[maxIndex];
-  Int_t rtIndex[maxIndex];
-
-  for(Int_t i = 0; i < maxIndex; i++) {
-    seedLabel[i] = kFALSE;
-    mcIndex[i] = -1;
-    rtIndex[i] = -1;
-  }
-  
-  // mark available seeds from TPC
-  Int_t nSeeds = 0;
-  Int_t nPrimarySeeds = 0;
-  
-  printf("marking found seeds from TPC\n"); 
-  TDirectory *savedir=gDirectory;  
-
-  TFile *in=TFile::Open("AliTPCBackTracks.root");
-  if (!in->IsOpen()) {
-    cerr<<"can't open file AliTPCBackTracks.root  !\n"; return;
-  }      
-  char   tname[100];
-  sprintf(tname,"seedsTPCtoTRD_%d",nEvent);
-  TTree *seedTree=(TTree*)in->Get(tname);
-  if (!seedTree) {
-     cerr<<"AliTRDtracker::PropagateBack(): ";
-     cerr<<"can't get a tree with seeds from TPC !\n";
-  }
-
-  AliTPCtrack *seed=new AliTPCtrack;
-  seedTree->SetBranchAddress("tracks",&seed);
-
-  Int_t n=(Int_t)seedTree->GetEntries();
-  for (Int_t i=0; i<n; i++) {
-     seedTree->GetEvent(i);
-     Int_t lbl = seed->GetLabel();
-     if(lbl < 0) { 
-       printf("negative seed label %d \n",lbl); 
-       continue;
-     }
-     if(lbl >= maxIndex) continue;
-     seedLabel[lbl] = kTRUE;
-     if(lbl < nPrimaries) nPrimarySeeds++;
-     nSeeds++;
-  }
-  delete seed;
-  delete seedTree;                                
-
-  printf("Found %d seeds from primaries among overall %d seeds \n",
-        nPrimarySeeds, nSeeds);
-
-  savedir->cd(); 
-  // done with marking TPC seeds
-
-
-  TFile *tf=TFile::Open("AliTRDtracks.root");
-
-  if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
-  TObjArray tarray(2000);
-
-  sprintf(tname,"TRDb_%d",nEvent);     
-  TTree *tracktree=(TTree*)tf->Get(tname);
-
-  TBranch *tbranch=tracktree->GetBranch("tracks");
-
-  Int_t nRecTracks = (Int_t) tracktree->GetEntries();
-  cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;
-
-  for (Int_t i=0; i<nRecTracks; i++) {
-    AliTRDtrack *iotrack=new AliTRDtrack();
-    tbranch->SetAddress(&iotrack);
-    tracktree->GetEvent(i);
-    tarray.AddLast(iotrack);
-    Int_t trackLabel = iotrack->GetLabel();
-
-    //    printf("rt with %d clusters and label %d \n",
-    //    iotrack->GetNumberOfClusters(), trackLabel);
-
-    if(trackLabel < 0) continue;
-    if(trackLabel >= maxIndex) continue;
-    rtIndex[trackLabel] = i;
-  }
-  tf->Close();                 
-
-  //  return;
-
-  // Load MC tracks 
-  TFile *mctf=TFile::Open("AliTRDmcTracks.root");
-  if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDmcTracks.root !\n"; return;}
-  TObjArray mctarray(2000);
-  TTree *mctracktree=(TTree*)mctf->Get("MCtracks");
-  TBranch *mctbranch=mctracktree->GetBranch("MCtracks");
-  Int_t nMCtracks = (Int_t) mctracktree->GetEntries();
-  cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
-  for (Int_t i=0; i<nMCtracks; i++) {
-    AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack;
-    mctbranch->SetAddress(&ioMCtrack);
-    mctracktree->GetEvent(i);
-    mctarray.AddLast(ioMCtrack);
-    Int_t mcLabel = ioMCtrack->GetTrackIndex();
-    if(mcLabel < 0) {printf("negative mc label detected!\n"); continue;}
-    if(mcLabel >= maxIndex) continue;
-    mcIndex[mcLabel] = i;
-  }
-  mctf->Close();                 
-
-
-  // Load clusters
-
-  TFile *geofile =TFile::Open("AliTRDclusters.root");   
-  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
-  Tracker->SetEventNumber(nEvent);
-
-  AliTRDgeometry *fGeom   = (AliTRDgeometry*) geofile->Get("TRDgeometry"); 
-  AliTRDparameter *fPar   = (AliTRDparameter*) geofile->Get("TRDparameter"); 
-  
-  Char_t *alifile = "AliTRDclusters.root";
-  TObjArray carray(2000);
-  TObjArray *ClustersArray = &carray;
-  Tracker->ReadClusters(ClustersArray,alifile);   
-
-
-  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
-  alifile = "galice.root";
-
-  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
-  if (!gafl) {
-    cout << "Open the ALIROOT-file " << alifile << endl;
-    gafl = new TFile(alifile);
-  }
-  else {
-    cout << alifile << " is already open" << endl;
-  }
-
-  // Get AliRun object from file or create it if not on file
-  gAlice = (AliRun*) gafl->Get("gAlice");
-  if (gAlice)
-    cout << "AliRun object found on file" << endl;
-  else
-    gAlice = new AliRun("gAlice","Alice test program");
-
-
-  // Define the objects
-  AliTRDv1       *TRD = (AliTRDv1*) gAlice->GetDetector("TRD");    
-
-  // Import the Trees for the event nEvent in the file
-  const Int_t nparticles = gAlice->GetEvent(nEvent);
-  if (nparticles <= 0) return;
-
-  TParticle *p;
-  Bool_t electrons[300000] = { kFALSE };
-
-  Bool_t mark_electrons = kFALSE;
-  if(mark_electrons) {
-    printf("mark electrons\n");
-
-    for(Int_t i = nPrimaries; i < nparticles; i++) {
-      p = gAlice->Particle(i);
-      if(p->GetMass() > 0.01) continue;
-      if(p->GetMass() < 0.00001) continue;
-      electrons[i] = kTRUE;
-    }
-  }
-
-  AliTRDcluster *cl = 0;
-  Int_t nw, label, index, ti[3], tbwc;
-  Int_t det, plane, ltb, gtb, gap, max_gap, sector, mc_sector, min_tb, max_tb;
-  Double_t Pt, Px, Py, Pz;
-  Double_t mcPt, mcPx, mcPy, mcPz;
-  Double_t x,y,z, mcX;
-  Int_t rtClusters, rtCorrect;
-
-  printf("\n");
-
-  AliTRDmcTrack *mct = 0;
-  AliTRDtrack *rt = 0;
-
-  Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
-  Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
-  Double_t dx = (Double_t) fPar->GetTimeBinSize();   
-
-  Int_t tbAmp = fPar->GetTimeBefore();
-  Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
-  Int_t tbDrift = fPar->GetTimeMax();
-  Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
-
-  tbDrift = TMath::Min(tbDrift,maxDrift);
-  tbAmp = TMath::Min(tbAmp,maxAmp);
-
-  const Int_t nPlanes = fGeom->Nplan();
-  const Int_t tbpp = Tracker->GetTimeBinsPerPlane();
-  const Int_t nTB = tbpp * nPlanes;
-  Int_t mask[nTB];
-
-  for(Int_t i=0; i < maxIndex; i++) {
-
-    rt = 0; mct = 0;
-    if(rtIndex[i] >= 0) rt = (AliTRDtrack *) tarray.UncheckedAt(rtIndex[i]);
-    if(mcIndex[i] >= 0) mct = (AliTRDmcTrack *) mctarray.UncheckedAt(mcIndex[i]);
-
-    if(!mct) continue;
-    label = mct->GetTrackIndex();
-    
-    //    if(TMath::Abs(mct->GetMass()-0.136) > 0.01) continue;
-
-    Int_t ncl = mct->GetNumberOfClusters();
-
-    // check how many time bins have a cluster
-
-    for(nw = 0; nw < nTB; nw++) mask[nw] = -1;
-
-    for(Int_t j = 0; j < ncl; j++) {
-      index = mct->GetClusterIndex(j);
-      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
-
-      for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
-
-      if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
-       printf("wrong track label: %d, %d, %d != %d \n",
-              ti[0], ti[1], ti[2], label);
-      } 
-      det=cl->GetDetector();
-      plane = fGeom->GetPlane(det); 
-      ltb = cl->GetLocalTimeBin();
-      gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-      mask[gtb] = index;
-    }
-
-
-
-
-    for(plane = 0; plane < nPlanes; plane++) {
-      for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
-       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-       if(mask[gtb] > -1) break;
-      }
-      if(ltb >= -tbAmp) break;
-    }  
-    if((plane == nPlanes) && (ltb == -tbAmp-1)) {
-      printf("warning: for track %d min tb is not found and set to %d!\n",
-            label, nTB-1);
-      min_tb = nTB-1;
-      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
-    }
-    else {
-      min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-    }
-
-    for(plane = nPlanes-1 ; plane>=0; plane--) {
-      for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
-       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-       if(mask[gtb] > -1) break;
-      }
-      if(ltb < tbDrift) break;
-    }  
-    if((plane == -1) && (ltb == tbDrift)) {
-      printf("warning: for track %d max tb is not found and set to 0!\n",label);
-      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
-      max_tb = 0;
-      mcX = Tracker->GetX(0,0,tbDrift-1);
-    }
-    else {
-      max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
-      mcX = Tracker->GetX(0,plane,ltb);
-      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb]);
-      det=cl->GetDetector();
-      sector = fGeom->GetSector(det);  
-      mc_sector = sector;
-    }
-
-
-    tbwc = 0;
-    max_gap = 0;
-    gap = -1;
-    
-    for(nw = min_tb; nw < max_tb+1; nw++) {
-      gap++;
-      if(mask[nw] > -1) {
-       tbwc++;
-       if(gap > max_gap) max_gap = gap;
-       gap = 0;
-      }
-    }
-
-    if(tbwc < ((Int_t) (nTB * Tracker->GetMinClustersInTrack()))) continue;
-    if(gap > Tracker->GetMaxGap()) continue; 
-    p = gAlice->Particle(label);
-    
-    printf("good track %d with min_tb %d, max_tb %d, gap %d\n",
-          label,min_tb,max_tb,gap);
-
-    hgood->Fill(p->Pt());
-
-    if(!rt) continue;
-
-    if(rt->GetLabel() != label) {
-      printf("mct vs rt index mismatch: %d != %d \n",
-            label, rt->GetLabel());
-      return;
-    }
-
-    if(!seedLabel[i]) continue;
-
-    hGoodAndSeed->Fill(p->Pt());
-
-    hXvsMCX->Fill(mcX,rt->GetX());
-    
-    rtClusters = rt->GetNumberOfClusters();
-
-    // find number of tb with correct cluster
-    rtCorrect = 0;
-    
-    for(Int_t j = 0; j < rtClusters; j++) {
-      index = rt->GetClusterIndex(j);
-      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
-
-      for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
-
-      if((ti[0] != label)&&(ti[1] != label)&&(ti[2] != label)) continue; 
-      rtCorrect++;
-    }
-    
-    Float_t  foundFraction = ((Float_t) rtCorrect) / (tbwc + 0.00001);
-    Float_t  correctFraction = ((Float_t) rtCorrect) / ((Float_t) rtClusters);
-
-    if(foundFraction > 1) printf("fraction = %f/%f \n",
-                                (Float_t) rtCorrect, (Float_t) tbwc);
-
-    
-    if(foundFraction <= 1) hFraction->Fill(foundFraction);
-    hCorrect->Fill(correctFraction);
-
-    hRTvsMC->Fill((Float_t) tbwc, (Float_t) rtClusters);
-    
-    if((foundFraction < 0.7) || (correctFraction < 0.7)) {
-      printf("not found track %d with FrctnFound %f and FrctnCorrect %f\n",
-            label, foundFraction, correctFraction);
-      continue;
-    }
-
-    hfound->Fill(p->Pt());
-
-    Pt = rt->Pt(); 
-    Double_t cc = TMath::Abs(rt->GetSigmaC2()); 
-    mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1);
-    rt->GetPxPyPz(Px,Py,Pz);      
-
-    printf("\n\ntrack %d \n", label);
-    printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz); 
-    printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz); 
-    
-    mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy);
-    if(mcPt < 0.0001) mcPt = 0.0001;
-    
-    Float_t lamg=TMath::ATan2(mcPz,mcPt);
-    Float_t lam =TMath::ATan2(Pz,Pt);
-    if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
-    else hpta->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc));
-    
-    Float_t phig=TMath::ATan2(mcPy,mcPx);
-    Float_t phi =TMath::ATan2(Py,Px);
-    
-
-    if(!(rt->PropagateTo(x))) continue;
-
-    
-    if((mcPt > Pt_min) && (mcPt < Pt_max)) {
-      hl->Fill((lam - lamg)*1000.);
-      hla->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2()));
-      if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
-      else hpt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.); 
-      hp->Fill((phi - phig)*1000.);
-      hy->Fill(rt->GetY() - y);
-      hz->Fill(rt->GetZ() - z);
-      hya->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2()));
-      hza->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2()));
-      hx->Fill((Float_t) x);
-    }  
-  }
-
-
-  gStyle->SetOptStat(0);
-  //  gStyle->SetOptStat(111110);
-  gStyle->SetOptFit(1); 
-
-
-  /*  
-  TCanvas *c1=new TCanvas("c1","",0,0,700,850);
-
-  //  gStyle->SetOptStat(0);
-  //  gStyle->SetOptStat(111110);
-  gStyle->SetOptFit(1); 
-    
-  TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
-  p1->cd(); p1->SetFillColor(10); p1->SetFrameFillColor(10);
-  hgood->Draw(); hGoodAndSeed->Draw("same"); c1->cd();  
-  
-  TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
-  p2->cd(); p2->SetFillColor(10); p2->SetFrameFillColor(10);
-  hFraction->SetYTitle("Counts");
-  hFraction->SetXTitle("Fraction");
-  hFraction->Draw(); c1->cd();      
-  
-  TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
-  p3->cd(); p3->SetFillColor(10); p3->SetFrameFillColor(10);
-  hfound->Draw(); c1->cd();   
-    
-  TPad *p4=new TPad("p4","",0.5,0,1.,0.3); p4->Draw();
-  p4->cd(); p4->SetFillColor(10); p4->SetFrameFillColor(10);
-  hCorrect->SetYTitle("Counts");
-  hCorrect->SetXTitle("Fraction");
-  hCorrect->Draw(); c1->cd();
-
-  TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
-  p5->SetFillColor(10); p5->SetFrameFillColor(10);
-  hgood->Sumw2(); hGoodAndSeed->Sumw2(); // hfake->Sumw2();
-  hSeedEff->Divide(hGoodAndSeed,hgood,1,1.,"b");
-  //  hf->Divide(hfake,hgood,1,1.,"b");
-  hSeedEff->SetMaximum(1.4);
-  hSeedEff->SetYTitle("TPC efficiency");
-  hSeedEff->SetXTitle("Pt (GeV/c)");
-  hSeedEff->Draw();          
-
-  TLine *line1 = new TLine(0,1.0,maxPt,1.0); line1->SetLineStyle(4);
-  line1->Draw("same");
-  TLine *line2 = new TLine(0,0.9,maxPt,0.9); line2->SetLineStyle(4);
-  line2->Draw("same");    
-  */  
-
-  TCanvas *c2=new TCanvas("c2","",0,0,700,850);
-    
-  TPad *p12=new TPad("p12","",0,0.3,.5,.6); p12->Draw();
-  p12->cd(); p12->SetFillColor(42); p12->SetFrameFillColor(10);
-  hp->SetFillColor(4);  hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c2->cd();  
-  
-  TPad *p22=new TPad("p22","",0.5,.3,1,.6); p22->Draw();
-  p22->cd(); p22->SetFillColor(42); p22->SetFrameFillColor(10);
-  hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c2->cd();      
-  
-  TPad *p32=new TPad("p32","",0,0,0.5,0.3); p32->Draw();
-  p32->cd(); p32->SetFillColor(42); p32->SetFrameFillColor(10);
-  hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c2->cd();   
-    
-  TPad *p42=new TPad("p42","",0.5,0,1.,0.3); p42->Draw();
-  p42->cd(); p42->SetFillColor(42); p42->SetFrameFillColor(10);
-  hgood->Draw(); hGoodAndSeed->Draw("same"); c2->cd();  
-
-  TPad *p52=new TPad("p52","",0,0.6,1,1); p52->Draw(); p52->cd();
-  p52->SetFillColor(41); p52->SetFrameFillColor(10);
-  hfound->Sumw2();
-  heff->Divide(hfound,hGoodAndSeed,1,1.,"b");
-  //  hf->Divide(hfake,hgood,1,1.,"b");
-  heff->SetMaximum(1.4);
-  heff->SetXTitle("Pt (GeV/c)");
-  heff->Draw();          
-
-  TLine *line12 = new TLine(0,1.0,maxPt,1.0); line12->SetLineStyle(4);
-  line12->Draw("same");
-  TLine *line22 = new TLine(0,0.9,maxPt,0.9); line22->SetLineStyle(4);
-  line22->Draw("same");    
-
-  c2->Print("matching.ps","ps"); 
-
-
-  /*  
-  TCanvas *cxyz = new TCanvas("cxyz","",50,50,750,900);
-  cxyz->Divide(2,2);
-  cxyz->cd(1); hx->Draw();
-  cxyz->cd(2); hy->Draw();
-  cxyz->cd(3); hz->Draw();
-  cxyz->cd(4); hXvsMCX->Draw();
-  */
-
-  TCanvas *cs = new TCanvas("cs","",0,0,700,850);
-  cs->Divide(2,3);
-
-  cs->cd(1); hy->Fit("gaus");
-  cs->cd(2); hz->Fit("gaus");
-  cs->cd(3); hpta->Fit("gaus");
-  cs->cd(4); hla->Fit("gaus");
-  cs->cd(5); hya->Fit("gaus");
-  cs->cd(6); hza->Fit("gaus");
-
-  cs->Print("resolution.ps","ps"); 
-
-  /*
-  TCanvas *cvs = new TCanvas("cvs","",0,0,700,850);
-  cvs->cd(); hRTvsMC->Draw();
-  */
-
-}
-
-
diff --git a/TRD/Macros/AliTRDclusterErrors.C b/TRD/Macros/AliTRDclusterErrors.C
deleted file mode 100644 (file)
index b63499d..0000000
+++ /dev/null
@@ -1,912 +0,0 @@
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <Riostream.h>
-
-#include "AliTRDtracker.h"
-#include "AliTRDclusterMI.h"
-#include "AliTRDhit.h"
-#include "AliTRDv1.h"
-#include "AliTRDgeometry.h"
-#include "AliTRDgeometryDetail.h"
-#include "AliTRDparameter.h"
-#include "AliTRDclusterCorrection.h"
-#include "alles.h"
-#include "TFile.h"
-#include "TStopwatch.h"
-#include "Rtypes.h"
-#include "TTree.h"
-
-#include "AliRunLoader.h"
-#include "AliStack.h"
-#include "TF1.h"
-#include "AliTrackReference.h"
-#endif    
-#include "AliTRDclusterErrors.h"
-
-
-AliTRDclusterCorrection * gCorrection;
-
-void ReadCorrection(){
-  TFile f("TRDcorrection.root");
-  gCorrection= (AliTRDclusterCorrection *)f.Get("TRDcorrection");
-  if (gCorrection==0){
-    printf("Correction not found");
-  }
-}
-
-
-class AliTRDExactPoint: public TObject {
-  public : 
-  AliTRDExactPoint();
-  Float_t fTX;    //x in rotated coordinate in the center of time bin
-  Float_t fTY;    //y 
-  Float_t fTZ;    //z
-  Float_t fTAY;   //angle y
-  Float_t fTAZ;   //angle z
-  Float_t fGx;
-  Float_t fGy;
-  Float_t fGz;
-  //
-  void SetReference(AliTrackReference *ref);
-  Float_t fTRefAngleY;
-  Float_t fRefPos[3];
-  Float_t fRefMom[3];
-  //
-  Int_t   fDetector;      // detector (chamber)
-  Int_t   fLocalTimeBin;  // local time bin
-  Int_t   fPlane;         // plane (layer)
-  Int_t   fSector;       // segment
-  Int_t   fPlaneMI;  
-  // 
-  Float_t fTQ;
-  Float_t fTPrim;
-  //  
-  ClassDef(AliTRDExactPoint,1)
-};
-
-class AliTRDCI: public TObject {
-  public :
-  AliTRDCI(){;}
-  virtual ~AliTRDCI(){;}
-  //
-  AliTRDclusterMI fCl;
-  AliTRDExactPoint fEp;
-  TParticle fP;
-  Char_t fStatus;
-  Int_t  fNClusters;
-  //
-  Float_t fDYtilt;
-  Float_t fTDistZ;  
-  //
-  Int_t   fNTracks;
-  Float_t fPt;
-  Float_t fCharge;
-  Bool_t  fIsPrim;
-  Float_t fCorrection;
-  void Update();
-  ClassDef(AliTRDCI,1)
-};
-
-class AliTRDClusterErrAnal: public TObject{
-public: 
-  AliTRDClusterErrAnal(Char_t *chloader  = "galice.root");
-  void SetIO(Int_t event);  
-  Int_t Analyze(Int_t trackmax);
-  void LoadClusters();
-  void MakeExactPoints(Int_t trackmax);
-  void SortReferences();
-  AliTrackReference * FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax=10.);
-public:
-  AliRunLoader * fRunLoader;
-  AliLoader * fTRDLoader;
-  AliTRDparameter *fParam;
-  AliTRDgeometry *fGeometry;
-  TTree * fHitTree;
-  TTree * fClusterTree;
-  TTree * fReferenceTree;
-  AliTRDv1 * fTRD;
-  //
-  TTree * fTreeA;
-  TFile * fFileA;
-  AliTRDtracker *fTracker;
-  AliStack *fStack;
-  TObjArray fClusters[12][100][18];  //first plane, second time bin
-  TObjArray fExactPoints;
-  TObjArray *fReferences;
-
-  ClassDef(AliTRDClusterErrAnal,1)
-};
-
-
-class AliTRDClusterErrDraw: public TObject{
-public:
-  TTree * fTree;
-  AliTRDclusterCorrection*   MakeCorrection(TTree * tree, Float_t offset);
-  static TH1F * ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin=10, Float_t ampmax=300);
-  static TH1F * ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min=-0.5, Float_t max=0.5);
-  static TH1F * ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min=-1., Float_t max=1.);
-  static void   AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle);
-  static TH1F*  CreateEffHisto(TH1F* hGen, TH1F* hRec);
-  static TH1F*  CreateResHisto(TH2F* hRes2, Bool_t draw = kTRUE, Bool_t drawBinFits = kTRUE, 
-                    Bool_t overflowBinFits = kFALSE);
-  ClassDef(AliTRDClusterErrDraw,1)
-};
-
-
-
-
-AliTRDExactPoint::AliTRDExactPoint()
-{
-  fTX=fTY=fTZ=fTAZ=fTAY=fGx=fGy=fGz=fTRefAngleY=0;
-  fRefPos[0]=fRefPos[1]=fRefPos[2]=fRefMom[0]=fRefMom[1]=fRefMom[2]=0;
-  fDetector=fLocalTimeBin=fPlane=fSector=fPlaneMI=0;
-  fTQ=fTPrim=0;
-}
-
-void AliTRDExactPoint::SetReference(AliTrackReference *ref){
-  fRefPos[0] = ref->X();
-  fRefPos[1] = ref->Y();
-  fRefPos[2] = ref->Z();
-  //
-  fRefMom[0] = ref->Px();
-  fRefMom[1] = ref->Py();
-  fRefMom[2] = ref->Pz();
-}
-
-
-void AliTRDCI::Update()
-{
-  //
-  //thanks to root
-  fPt = fP.Pt();
-  fCharge = fP.GetPDG()->Charge();
-  fIsPrim = (fP.GetMother(0)>0)? kFALSE :kTRUE;
-  fCorrection = gCorrection->GetCorrection(fEp.fPlane,fCl.fTimeBin,fEp.fTAY);
-}
-
-
-/*
-//example seesion
-
-.L AliGenInfo.C+
-.L AliTRDclusterErrors.C+
-gCorrection  = AliTRDclusterCorrection::GetCorrection();
-AliTRDClusterErrAnal ana;
-ana.Analyze(10000000)
-
-
-
-.L AliGenInfo.C+
-.L AliTRDclusterErrors.C+
-TFile f("trdclusteranal.root");
-TTree* tree = (TTree*)f.Get("trdcl");
-AliComparisonDraw comp;
-comp->fTree = tree;
-
-
-tree->SetAlias("shapef","(1.-(0.8+0.06*(6-fEp.fPlane))*(fCl.fSigmaY2/(0.17+0.027*abs(fEp.fTAY))))");
-tree->SetAlias("shapes","0.08+0.3/sqrt(fCl.fQ)");
-tree->SetAlias("sfactor","shapef/shapes");
-
-tree->SetAlias("shapen","(fCl.fNPads-2.7-0.9*abs(fEp.fTAY))");
-
-tree->SetAlias("gshape","sfactor>-2&&fCl.fNPads<6&&shapen<1");
-
-
-tree->SetAlias("dy"    , "fEp.fTY-fCl.fY-fDYtilt");
-tree->SetAlias("angle","abs(fEp.fTAY)");
-TCut cbase("cbase","(abs(fP.fPdgCode)!=11||fPt>0.2)&&fPt>0.1&&angle<2");
-
-tree->SetAlias("erry0","(0.028+0.07*angle)");
-tree->SetAlias("erry1","erry0*(0.9+15./fCl.fQ)");
-tree->SetAlias("erry2","erry1*(0.8+0.5*max(-sfactor,0))"); 
-
-
-
-
-TH1F his("resy","resy",100,-0.2,0.2);
-comp->fTree->Draw("dy:0.028*fEp.fTAY","fStatus==0 && abs(dy)<1.0&&fNTracks<2&&angle<2&&gshape","")
-comp->DrawXY("sqrt(fCl.fQ)","dy","fStatus==0"+cbase,"gshape",10,0,20,-0.7,0.7);
-
-comp->DrawXY("angle","dy/erry1","fStatus==0"+cbase,"gshape",10,0,2,-5.7,5.7);
-comp->DrawXY("sqrt(cl->fQ)","dy/err1","fStatus==0"+cbase,"gshape",10,2,20,-5.7,5.7);
-
-
-
-AliTRDClusterErrDraw::ResDyVsAmp(tree,"abs(dy)<0.4&&abs(fEp.fTAY)<0.2&&fNTracks<2&&fPt>0.1"+cbase,-0.03);
-AliTRDClusterErrDraw::ResDyVsRelPos(tree,"abs(dy)<0.4&&abs(fEp.fTAY)<0.2&&fNTracks<2&&fPt>0.2"+cbase,-0.03);
-AliTRDClusterErrDraw::ResDyVsAngleY(tree,"abs(fEp.fTY-fCl.fY+fDYtilt)<0.4",0.);
-
-AliTRDclusterCorrection * cor = AliTRDClusterErrDraw::MakeCorrection(tree,0.134);
-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");
-
-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") 
-
-
-*/
-
-
-ClassImp(AliTRDCI)
-ClassImp(AliTRDExactPoint)
-ClassImp(AliTRDClusterErrAnal)
-ClassImp(AliTRDClusterErrDraw)
-
-
-
-AliTRDClusterErrAnal::AliTRDClusterErrAnal(Char_t *chloader )
-{
-  //
-  //SET Input loaders
-  if (gAlice){
-    delete gAlice->GetRunLoader();
-    delete gAlice;
-    gAlice = 0x0;
-  }    
-  fRunLoader = AliRunLoader::Open(chloader);
-  if (fRunLoader == 0x0){
-    cerr<<"Can not open session"<<endl;
-    return;
-  }
-  fTRDLoader = fRunLoader->GetLoader("TRDLoader");  
-  if (fTRDLoader == 0x0){
-    cerr<<"Can not get TRD Loader"<<endl;
-    return ;
-  }   
-  if (fRunLoader->LoadgAlice()){
-    cerr<<"Error occured while l"<<endl;
-    return;
-  }
-  fRunLoader->CdGAFile();
-  fTracker = new AliTRDtracker(gFile);
-  fParam = (AliTRDparameter*) gFile->Get("TRDparameter");
-  fGeometry = new AliTRDgeometryDetail();   
-  fTRD      = (AliTRDv1*) gAlice->GetDetector("TRD");
-  //
-  AliTRDCI * clinfo = new AliTRDCI();
-  fFileA  = new TFile("trdclusteranal.root","recreate");
-  fFileA->cd();
-  fTreeA  = new TTree("trdcl","trdcl");
-  fTreeA->Branch("trdcl","AliTRDCI",&clinfo);
-  delete clinfo;
-}
-
-void AliTRDClusterErrAnal::SetIO(Int_t event)
-{
-  //
-  //set input output for given event
-  fRunLoader->SetEventNumber(event);
-  fRunLoader->LoadHeader();
-  fRunLoader->LoadKinematics();
-  fRunLoader->LoadTrackRefs();
-  fTRDLoader->LoadHits();
-  fTRDLoader->LoadRecPoints("read");
-  //
-  fStack = fRunLoader->Stack();
-  fHitTree = fTRDLoader->TreeH();
-  fClusterTree = fTRDLoader->TreeR();
-  fReferenceTree = fRunLoader->TreeTR();
-  fTracker->LoadClusters(fClusterTree);
-  //
-}
-
-void AliTRDClusterErrAnal::LoadClusters()
-{
-  //
-  //
-  // Load clusters  
-  TObjArray *ClusterArray = new TObjArray(400);   
-  TObjArray carray(2000);
-  TBranch *branch=fClusterTree->GetBranch("TRDcluster");
-  if (!branch) {
-    Error("ReadClusters","Can't get the branch !");
-    return;
-  }
-  Int_t over5 =0;
-  Int_t over10=0;
-
-  branch->SetAddress(&ClusterArray);
-  Int_t nentries = (Int_t)fClusterTree->GetEntries();
-  for (Int_t i=0;i<nentries;i++){
-    fClusterTree->GetEvent(i);
-    Int_t nCluster = ClusterArray->GetEntriesFast();
-    for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
-      AliTRDcluster* c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
-      carray.AddLast(c);
-      ClusterArray->RemoveAt(iCluster);
-      if (c->GetQ()>5)  over5++;
-      if (c->GetQ()>10) over10++;
-    }
-  }
-  Int_t nClusters = carray.GetEntriesFast();
-  printf("Total number of clusters %d\t%d\t%d\n", nClusters,over5,over10);
-  //
-  //
-  //SORT clusters
-  //
-  Int_t all=0;  
-  for (Int_t i=0;i<nClusters;i++){
-    AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);  
-    Int_t plane = fGeometry->GetPlane(cl->GetDetector());
-    if (plane>=12) continue;
-    Int_t time = cl->GetLocalTimeBin(); 
-    if (time>=100) continue;
-    Int_t sector = fGeometry->GetSector(cl->GetDetector());
-    if (sector>=18){
-      printf("problem1\n");
-      continue;
-    }    
-    fClusters[plane][time][sector].AddLast(cl);
-    all++;
-  }
-}
-
-void AliTRDClusterErrAnal::SortReferences()
-{
-  //
-  //
-  //
-  printf("Sorting references\n");
-  fReferences = new TObjArray;
-  Int_t ntracks = fStack->GetNtrack();
-  fReferences->Expand(ntracks);
-  Int_t nentries = (Int_t)fReferenceTree->GetEntries();
-  TClonesArray * arr = new TClonesArray("AliTrackReference");
-  TBranch * br = fReferenceTree->GetBranch("TRD");
-  br->SetAddress(&arr);
-  //
-  TClonesArray *labarr=0;
-  Int_t nreferences=0;
-  Int_t nreftracks=0;
-  for (Int_t iprim=0;iprim<nentries;iprim++){
-    if (br->GetEntry(iprim)){
-      for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
-       AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
-       if (!ref) continue;
-       Int_t lab = ref->GetTrack();
-       if ( (lab<0) || (lab>ntracks)) continue;
-       //
-       if (fReferences->At(lab)==0) {
-         labarr = new TClonesArray("AliTrackReference"); 
-         fReferences->AddAt(labarr,lab);
-         nreftracks++;
-       }
-       TClonesArray &larr = *labarr;
-       new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
-       nreferences++;
-      }
-    }
-  }
-  printf("Total number of references = \t%d\n", nreferences);
-  printf("Total number of tracks with references = \t%d\n", nreftracks);
-  printf("End - Sorting references\n");
-  
-}
-
-AliTrackReference * AliTRDClusterErrAnal::FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax)
-{
-  //
-  //
-  //
-  if (fReferences->At(lab)==0) return 0;
-  AliTrackReference *nearest=0;
-  TClonesArray * arr = (TClonesArray *)fReferences->At(lab);
-  for (Int_t iref =0;iref<arr->GetEntriesFast();iref++){
-    AliTrackReference * ref = ( AliTrackReference *)arr->UncheckedAt(iref);
-    if (!ref) continue;
-    Float_t delta = (pos[0]-ref->X())*(pos[0]-ref->X());
-    delta += (pos[1]-ref->Y())*(pos[1]-ref->Y());
-    delta += (pos[2]-ref->Z())*(pos[2]-ref->Z());
-    delta = TMath::Sqrt(delta);
-    if (delta<dmax){
-      dmax=delta;
-      nearest = ref;
-    }
-  }
-  return nearest;
-}
-
-void AliTRDClusterErrAnal::MakeExactPoints(Int_t trackmax)
-{
-  //
-  //make exact points:-)       
-  //
-  // 
-
-  fExactPoints.Delete();
-  fExactPoints.Expand(fStack->GetNtrack());
-  //
-  Double_t fSum=0;
-  Double_t fSumQ =0;
-  Double_t fSumX=0;
-  Double_t fSumX2=0;
-  Double_t fSumXY=0;
-  Double_t fSumXZ=0;
-  Double_t fSumY=0;
-  Double_t fSumZ=0;
-  //
-  Int_t entries = Int_t(fHitTree->GetEntries());
-  printf("Number of primary  entries\t%d\n",entries);
-  entries = TMath::Min(trackmax,entries);
-  Int_t nallpoints = 0;
-
-  Int_t nalltracks =0;
-  Int_t pointspertrack =0;
-
-  for (Int_t entry=0;entry<entries; entry++){
-    gAlice->ResetHits();
-    fHitTree->GetEvent(entry);
-    Int_t lastlabel    = -1;
-    Int_t lastdetector = -1;
-    Int_t lasttimebin   = -1;
-    Float_t lastpos[3];
-    //
-    for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); hit; 
-       hit = (AliTRDhit *) fTRD->NextHit()) {
-      //
-      Int_t label    = hit->Track();
-      TParticle * particle = fStack->Particle(label);
-      if (!particle) continue;
-      if (particle->Pt()<0.05) continue;
-      Int_t detector = hit->GetDetector();
-      Int_t plane    = fGeometry->GetPlane(detector);
-      //
-      //
-      if (hit->GetCharge()==0) continue;
-      Float_t pos[3] = {hit->X(),hit->Y(),hit->Z()};
-      Int_t indexes[3];
-      fGeometry->Global2Detector(pos,indexes,fParam);
-      //
-      Float_t rot[3];
-      fGeometry->Rotate(detector,pos,rot);
-      //rot[0] *=-1;
-      //  rot[1] *=-1;
-      //
-      //
-      Float_t  time0    = fParam->GetTime0(plane);
-      Int_t    timebin  = Int_t(TMath::Nint(((time0 - rot[0])/fParam->GetTimeBinSize())+ fParam->GetTimeBefore())+0.1); 
-      if (timebin<0) continue;
-      //
-      //
-      if (label!=lastlabel || detector != lastdetector || lasttimebin !=timebin){
-       //
-       if (label!=lastlabel){
-         fExactPoints.AddAt(new TClonesArray("AliTRDExactPoint",0),label);
-         //printf("new particle\t%d\n",label);
-         nalltracks++;
-         //      printf("particle\t%d- hits\t%d\n",lastlabel, pointspertrack);
-         pointspertrack=0;
-       }
-       
-       if ( (fSum>1) && lasttimebin>=0 && lasttimebin<fParam->GetTimeMax() ){
-         //if we have enough info for given layer time bin - store it
-         AliTrackReference * ref = FindNearestReference(lastlabel,lastpos,4.);
-         Float_t rotmom[3];
-         Float_t rotpos[3];
-         Float_t refangle=0;
-         if (ref){
-           Float_t mom[3] = {ref->Px(),ref->Py(),ref->Pz()};
-           Float_t refpos[3] = {ref->X(),ref->Y(),ref->Z()};
-           fGeometry->Rotate(detector,mom,rotmom);
-           fGeometry->Rotate(detector,refpos,rotpos);
-           refangle = rotmom[1]/rotmom[0];
-
-         }
-
-         Double_t ay,by,az,bz;
-         Double_t det = fSum*fSumX2-fSumX*fSumX;
-         if (TMath::Abs(det)> 0.000000000000001) { 
-           by = (fSum*fSumXY-fSumX*fSumY)/det;
-           ay = (fSumX2*fSumY-fSumX*fSumXY)/det;
-           
-         }else{
-           ay =fSumXY/fSumX;
-           by =0;         
-         }
-         if (TMath::Abs(det)> 0.000000000000001) { 
-           bz = (fSum*fSumXZ-fSumX*fSumZ)/det;
-           az = (fSumX2*fSumZ-fSumX*fSumXZ)/det;         
-         }else{
-           az =fSumXZ/fSumX;
-           bz =0;         
-         }
-         //
-         Float_t lastplane = fGeometry->GetPlane(lastdetector);
-         Float_t time0    = fParam->GetTime0(lastplane);
-         Float_t xcenter0 = time0 - (lasttimebin - fParam->GetTimeBefore()+0.5)*fParam->GetTimeBinSize();
-         Float_t xcenter = fTracker->GetX(0,lastplane,lasttimebin);
-         if (TMath::Abs(xcenter-xcenter0)>0.001){
-           printf("problem");
-         }
-         
-         Float_t ty = ay + by * xcenter;
-         Float_t tz = az + bz * xcenter;
-         //
-
-         TClonesArray * arr = (TClonesArray *) fExactPoints.At(label);
-         TClonesArray & larr= *arr;
-         Int_t arrent = arr->GetEntriesFast();
-         AliTRDExactPoint * point = new (larr[arrent]) AliTRDExactPoint;
-         nallpoints++;
-
-         if (ref){
-           point->SetReference(ref);
-           point->fTRefAngleY = rotmom[1]/rotmom[0];
-         }
-         point->fTX  = xcenter;
-         point->fTY  = ty;
-         point->fTZ  = tz;
-         point->fTAY = by;
-         point->fTAZ = bz;
-         //
-         point->fGx  = lastpos[0];
-         point->fGy  = lastpos[1];
-         point->fGz  = lastpos[2];
-
-         //
-         point->fDetector     = lastdetector;
-         point->fLocalTimeBin = lasttimebin;
-         point->fPlane        = fGeometry->GetPlane(lastdetector); 
-         point->fSector      = fGeometry->GetSector(lastdetector); 
-         point->fPlaneMI      = indexes[0]; 
-         //
-         point->fTPrim = fSum;
-         point->fTQ    = fSumQ;            
-         //
-       }
-       lastdetector = detector;
-       lastlabel    = label;
-       lasttimebin  = timebin;
-       fSum=fSumQ=fSumX=fSumX2=fSumXY=fSumXZ=fSumY=fSumZ=0.;
-      }
-      //
-      lastpos[0] = hit->X();
-      lastpos[1] = hit->Y();
-      lastpos[2] = hit->Z();
-      fSum++;
-      fSumQ  +=hit->GetCharge();
-      fSumX  +=rot[0];
-      fSumX2 +=rot[0]*rot[0];
-      fSumXY +=rot[0]*rot[1];      
-      fSumXZ +=rot[0]*rot[2];
-      fSumY  +=rot[1];      
-      fSumZ  +=rot[2];     
-      pointspertrack++;
-    }
-  }
-  //
-  printf("Found %d exact points\n",nallpoints); 
-}
-
-
-
-
-
-
-Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {  
-
-  //
-  // comparison works with both cluster types MI and old also
-  //dummy cluster to be fill if not cluster info
-  AliTRDclusterMI clmi;  
-  TClass * classmi =  clmi.IsA();
-  //
-  //SetOutput
-  AliTRDCI * clinfo = new AliTRDCI();
-  TBranch * clbr = fTreeA->GetBranch("trdcl");
-  clbr->SetAddress(&clinfo);
-
-  SetIO(0);
-  SortReferences();
-  MakeExactPoints(trackmax);
-  LoadClusters();
-  //
-  trackmax =  fStack->GetNtrack();
-  //
-  // Get the number of entries in the hit tree
-  // (Number of primary particles creating a hit somewhere)
-  Int_t nTrack = (Int_t)fExactPoints.GetEntries();
-  printf("Found %d charged in TRD in first %d particles", nTrack, trackmax);
-  //
-
-  for (Int_t itrack = 0; itrack<trackmax; itrack++){
-    TClonesArray *arrpoints = (TClonesArray*)fExactPoints.At(itrack);
-    
-    if (!arrpoints) continue;
-    //printf("new particle\t%d\n",itrack);
-    TParticle * particle = fStack->Particle(itrack);
-    if (!particle) continue;
-    //printf("founded in kine tree \t%d\n",itrack);
-    Int_t npoints = arrpoints->GetEntriesFast();
-    if (npoints<10) continue;
-    //printf("have enough points \t%d\t%d\n",itrack,npoints);
-
-    for (Int_t ipoint=0;ipoint<npoints;ipoint++){
-      AliTRDExactPoint * point = (AliTRDExactPoint*)arrpoints->UncheckedAt(ipoint);
-      if (!point) continue;
-      //
-      Int_t sec = fGeometry->GetSector(point->fDetector);
-      if (sec>18){
-       printf("problem2\n");
-      }
-      TObjArray & cllocal = fClusters[point->fPlane][point->fLocalTimeBin][sec];       
-      Int_t nclusters = cllocal.GetEntriesFast();
-      Float_t maxdist = 10;
-      AliTRDcluster * nearestcluster =0;
-      clinfo->fNClusters=0;
-      //find nearest cluster to hit with given label
-      for (Int_t icluster =0; icluster<nclusters; icluster++){
-       AliTRDcluster * cluster = (AliTRDcluster*)cllocal.UncheckedAt(icluster);
-       if (!cluster) continue;
-       if ( (cluster->GetLabel(0)!= itrack) &&  (cluster->GetLabel(1)!= itrack)&&(cluster->GetLabel(2)!= itrack))
-         continue;
-       Float_t dist = TMath::Abs(cluster->GetY()-point->fTY);
-       if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5 || dist>3.) continue; 
-       clinfo->fNClusters++;
-       if (dist<maxdist){
-         maxdist = dist;
-         nearestcluster = cluster;
-       }
-      }      
-      //
-      clinfo->fEp  = *point;
-      clinfo->fP   = *particle;
-      if (!nearestcluster) {
-       clinfo->fStatus=1;
-       clinfo->fCl = clmi;
-      }
-      else{
-       clinfo->fStatus=0;
-       if (nearestcluster->IsA()==classmi){
-         clinfo->fCl    =*((AliTRDclusterMI*)nearestcluster);    
-       }
-       else{
-         clinfo->fCl    = *nearestcluster;
-       }
-       //     
-       Float_t dz      = clinfo->fCl.GetZ()-point->fTZ;
-       Double_t h01    = sin(TMath::Pi() / 180.0 * fParam->GetTiltingAngle());
-       clinfo->fTDistZ = dz;
-       clinfo->fDYtilt = h01*dz*((point->fPlane%2)*2.-1.); 
-       //
-       clinfo->fNTracks =1;
-       if (nearestcluster->GetLabel(1)>=0)  clinfo->fNTracks++;
-       if (nearestcluster->GetLabel(2)>=0)  clinfo->fNTracks++;  
-       clinfo->Update();
-      }
-      //
-      fTreeA->Fill();
-    }
-  }
-  
-  
-  fFileA->cd();
-  fTreeA->Write();
-  fFileA->Close();
-  return 0;
-}
-
-AliTRDclusterCorrection*   AliTRDClusterErrDraw::MakeCorrection(TTree * tree, Float_t offset)
-{
-  //
-  //
-  // make corrections
-  AliTRDclusterCorrection * cor = new AliTRDclusterCorrection;
-  cor->SetOffsetAngle(offset); 
-  for (Int_t iplane=0;iplane<6;iplane++)
-    for (Int_t itime=0;itime<15;itime++)
-      for (Int_t iangle=0; iangle<20;iangle++){
-       Float_t angle = cor->GetAngle(iangle);
-       TH1F delta("delta","delta",30,-0.3,0.3);
-       char selection[100]="fStatus==0&&fNTracks<2";
-       char selectionall[1000];
-       sprintf(selectionall,"%s&&abs(fEp.fTAY-%f)<0.2&&fEp.fPlane==%d&&fCl.fTimeBin==%d&&fCl.fQ>20",
-               selection,angle,iplane,itime);
-       printf("\n%s",selectionall);
-       tree->Draw("fEp.fTY-fCl.fY+fDYtilt>>delta",selectionall);
-       gPad->Update();
-       printf("\nplane\t%d\ttime%d\tangle%f",iplane,itime,angle);
-       printf("\tentries%f\tmean\t%f\t%f",delta.GetEntries(),delta.GetMean(),delta.GetRMS());  
-       cor->SetCorrection(iplane,itime,angle,delta.GetMean(),delta.GetRMS());
-      }
-  TFile * f = new TFile("TRDcorrection.root","new");
-  if (!f) f = new TFile("TRDcorrection.root","recreate");
-  f->cd();
-  cor->Write("TRDcorrection");
-  f->Close();
-  return cor; 
-}
-
-TH1F * AliTRDClusterErrDraw::ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin, Float_t ampmax)
-{
-  //
-  //
-  TH2F hisdy("resy","resy",10,ampmin,ampmax,30,-0.3,0.3);
-  char expression[1000];
-  sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fQ>>resy",t0);
-  char selectionall[1000];
-  sprintf(selectionall,"fStatus==0&&%s",selection);
-  printf("%s\n",expression);
-  printf("%s\n",selectionall);
-  tree->Draw(expression,selectionall);
-  return CreateResHisto(&hisdy);
-}
-
-
-TH1F * AliTRDClusterErrDraw::ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
-{
-  //
-  //
-  min *=128;
-  max *=128;
-  TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
-  char expression[1000];
-  sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fRelPos>>resy",t0);
-  char selectionall[1000];
-  sprintf(selectionall,"fStatus==0&&%s",selection);
-  printf("%s\n",expression);
-  printf("%s\n",selectionall);
-  tree->Draw(expression,selectionall);
-  return CreateResHisto(&hisdy);
-}
-
-
-TH1F * AliTRDClusterErrDraw::ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
-{
-  //
-  //
-  TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
-
-  char expression[1000];
-  sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%f*fEp.fTAY:fEp.fTAY>>resy",t0);
-  char selectionall[1000];
-  sprintf(selectionall,"fStatus==0&&%s",selection);
-
-  tree->Draw(expression,selectionall);
-  return CreateResHisto(&hisdy);
-}
-
-void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
-{
-  histo->GetXaxis()->SetTitle(xAxisTitle);
-  histo->GetYaxis()->SetTitle(yAxisTitle);
-}
-
-
-TH1F* AliTRDClusterErrDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
-{
-  Int_t nBins = hGen->GetNbinsX();
-  TH1F* hEff = (TH1F*) hGen->Clone("hEff");
-  hEff->SetTitle("");
-  hEff->SetStats(kFALSE);
-  hEff->SetMinimum(0.);
-  hEff->SetMaximum(110.);
-
-  for (Int_t iBin = 0; iBin <= nBins; iBin++) {
-    Double_t nGen = hGen->GetBinContent(iBin);
-    Double_t nRec = hRec->GetBinContent(iBin);
-    if (nGen > 0) {
-      Double_t eff = nRec/nGen;
-      hEff->SetBinContent(iBin, 100. * eff);
-      Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
-      //      if (error == 0) error = sqrt(0.1/nGen);
-      //
-      if (error == 0) error = 0.0001;
-      hEff->SetBinError(iBin, 100. * error);
-    } else {
-      hEff->SetBinContent(iBin, 100. * 0.5);
-      hEff->SetBinError(iBin, 100. * 0.5);
-    }
-  }
-
-  return hEff;
-}
-
-
-
-TH1F* AliTRDClusterErrDraw::CreateResHisto(TH2F* hRes2, Bool_t draw,  Bool_t drawBinFits, 
-                    Bool_t overflowBinFits)
-{
-  TVirtualPad* currentPad = gPad;
-  TAxis* axis = hRes2->GetXaxis();
-  Int_t nBins = axis->GetNbins();
-  TH1F* hRes, *hMean;
-  if (axis->GetXbins()->GetSize()){
-    hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
-    hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
-  }
-  else{
-    hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
-    hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
-
-  }
-  hRes->SetStats(false);
-  hRes->SetOption("E");
-  hRes->SetMinimum(0.);
-  //
-  hMean->SetStats(false);
-  hMean->SetOption("E");
-  // create the fit function
-  //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
-  TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
-  
-  fitFunc->SetLineWidth(2);
-  fitFunc->SetFillStyle(0);
-  // create canvas for fits
-  TCanvas* canBinFits = NULL;
-  Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
-  Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
-  Int_t ny = (nPads-1) / nx + 1;
-  if (drawBinFits) {
-    canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
-    if (canBinFits) delete canBinFits;
-    canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
-    canBinFits->Divide(nx, ny);
-  }
-
-  // loop over x bins and fit projection
-  Int_t dBin = ((overflowBinFits) ? 1 : 0);
-  for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
-    if (drawBinFits) canBinFits->cd(bin + dBin);
-    TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
-    //    
-    if (hBin->GetEntries() > 10) {
-      fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS(),0.02*hBin->GetMaximum());
-      hBin->Fit(fitFunc,"s");
-      Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
-
-      if (sigma > 0.){
-       hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
-       hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
-      }
-      else{
-       hRes->SetBinContent(bin, 0.);
-       hMean->SetBinContent(bin,0);
-      }
-      hRes->SetBinError(bin, fitFunc->GetParError(2));
-      hMean->SetBinError(bin, fitFunc->GetParError(1));
-      
-      //
-      //
-
-    } else {
-      hRes->SetBinContent(bin, 0.);
-      hRes->SetBinError(bin, 0.);
-      hMean->SetBinContent(bin, 0.);
-      hMean->SetBinError(bin, 0.);
-    }
-    
-
-    if (drawBinFits) {
-      char name[256];
-      if (bin == 0) {
-       sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
-      } else if (bin == nBins+1) {
-       sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
-      } else {
-       sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
-               axis->GetTitle(), axis->GetBinUpEdge(bin));
-      }
-      canBinFits->cd(bin + dBin);
-      hBin->SetTitle(name);
-      hBin->SetStats(kTRUE);
-      hBin->DrawCopy("E");
-      canBinFits->Update();
-      canBinFits->Modified();
-      canBinFits->Update();
-    }
-    
-    delete hBin;
-  }
-
-  delete fitFunc;
-  currentPad->cd();
-  if (draw){
-    currentPad->Divide(1,2);
-    currentPad->cd(1);
-    hRes->Draw();
-    currentPad->cd(2);
-    hMean->Draw();
-  }
-  
-  return hRes;
-}
diff --git a/TRD/Macros/AliTRDsaveMCtracks.C b/TRD/Macros/AliTRDsaveMCtracks.C
deleted file mode 100644 (file)
index 281405a..0000000
+++ /dev/null
@@ -1,314 +0,0 @@
-#ifndef __CINT__
-  #include <iostream.h>
-  #include "AliTRDtracker.h"
-
-  #include "AliTRDcluster.h" 
-  #include "AliTRDhit.h" 
-  #include "AliTRDv1.h"
-  #include "AliTRDgeometry.h"    
-  #include "AliTRDparameter.h"    
-  #include "alles.h"  
-  #include "AliTRDmcTrack.h"
-
-  #include "AliTRDtrack.h"
-  
-  #include "TFile.h"
-  #include "TParticle.h"
-  #include "TStopwatch.h"
-#endif
-
-Int_t Find(Int_t prtcl, TObjArray *mctarray) {
-
-// Returns index of the mct which corresponds to particle <prtcl> 
-
-  Int_t b=0, e=mctarray->GetEntriesFast(), m=(b+e)/2;
-  AliTRDmcTrack *mct = 0; 
-  while ((b<e)&&(e!=0)) {
-    m=(b+e)/2;
-    mct = (AliTRDmcTrack*) mctarray->UncheckedAt(m);
-    if (prtcl > mct->GetTrackIndex()) b=m+1;
-    else e=m;
-  }
-
-  if(mct->GetTrackIndex() == prtcl) return m;
-
-  if((m+1) < mctarray->GetEntriesFast()) {
-    mct = (AliTRDmcTrack*) mctarray->UncheckedAt(m+1);
-    if(mct->GetTrackIndex() == prtcl) return m+1;
-  }
-
-  return -1;
-}                    
-
-void AliTRDsaveMCtracks() {
-
-  TObjArray mctracks(2000);
-  TObjArray *TRDmcTracks = &mctracks; 
-
-  Char_t *alifile = "galice.root";
-  Int_t   nEvent  = 0; 
-  Float_t low_pt_cut = 0.05;
-  Float_t low_p_cut = 0.02;
-  Int_t min_no_of_clusters = 10;
-
-  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
-  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
-  if (!gafl) {
-    cout << "Open the ALIROOT-file " << alifile << endl;
-    gafl = new TFile(alifile);
-  }
-  else {
-    cout << alifile << " is already open" << endl;
-  }
-
-  // Get AliRun object from file or create it if not on file
-  gAlice = (AliRun*) gafl->Get("gAlice");
-  if (gAlice)
-    cout << "AliRun object found on file" << endl;
-  else
-    gAlice = new AliRun("gAlice","Alice test program");
-
-  AliTRDv1       *fTRD     = (AliTRDv1*) gAlice->GetDetector("TRD");   
-  
-  // Import the Trees for the event nEvent in the file
-  const Int_t nparticles = gAlice->GetEvent(nEvent);
-
-  printf("found %d particles in event %d \n", nparticles, nEvent);
-    
-  // Create TRDmcTracks for each particle
-  Int_t label = -1, charge = 0;
-  Bool_t primary;
-  Float_t mass;
-  Int_t pdg_code;
-
-  printf("\n");
-
-  for (Int_t ii=0; ii<nparticles; ii++) {
-
-    printf("\r particle %d out of %d",ii,nparticles);
-
-    TParticle *p = gAlice->Particle(ii);
-    if(p->P() < low_p_cut) continue;
-
-    primary = kTRUE;
-    if (p->GetFirstMother()>=0) primary = kFALSE;
-    
-    pdg_code = (Int_t) p->GetPdgCode();
-
-    if ((pdg_code == 10010020) ||
-       (pdg_code == 10010030) || 
-       (pdg_code == 50000050) || 
-       (pdg_code == 50000051) || 
-       (pdg_code == 10020040)) {
-
-      mass = 0.;
-      charge = 0;
-    }
-    else {
-      TParticlePDG *pdg = p->GetPDG();
-      charge = (Int_t) pdg->Charge(); 
-      mass=pdg->Mass();
-    }
-    if(charge == 0) continue;
-
-    AliTRDmcTrack *t = new AliTRDmcTrack(ii, primary, mass, charge, pdg_code);
-    TRDmcTracks->AddLast(t);
-  }
-  printf("\r\n");
-
-  // Loop through TRD clusters and assign indexes to MC tracks
-  
-  TFile *geofile =TFile::Open("AliTRDclusters.root");   
-  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
-  Tracker->SetEventNumber(nEvent);
-
-  AliTRDgeometry *fGeo   = (AliTRDgeometry*) geofile->Get("TRDgeometry"); 
-  AliTRDparameter *fPar   = (AliTRDparameter*) geofile->Get("TRDparameter"); 
-
-
-  Char_t *clusterfile = "AliTRDclusters.root";
-  TObjArray carray(2000);
-  TObjArray *ClustersArray = &carray;  
-  Tracker->ReadClusters(ClustersArray,clusterfile);
-
-  Int_t nClusters = carray.GetEntriesFast();
-
-  Int_t track, index;
-  AliTRDmcTrack *tpr = NULL;
-
-  for (Int_t i = 0; i < nClusters; i++) {
-
-    printf("\r assigning cluster %d out of %d", i, nClusters);
-    AliTRDcluster *cl = (AliTRDcluster *) ClustersArray->UncheckedAt(i);
-
-    for(Int_t j=0; j<3; j++) {
-      track = cl->GetLabel(j);
-      if(track < 0) continue;
-      if(track >= nparticles) 
-       printf("Track index %d is larger than total number of particles %d\n"
-               ,track,nparticles);
-      else {
-       index = Find(track, TRDmcTracks);
-       tpr = 0;
-       if(index > 0) tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index);
-       
-       if(tpr) {
-         label = tpr->GetTrackIndex();
-         if(label != track) printf("Got label %d while expected %d !\n",
-                                   label, track);
-         TParticle *p = gAlice->Particle(track);
-
-         if (p->Pt() > low_pt_cut) tpr->Update(i);       
-       }         
-      }
-    }
-  }
-  
-  // Loop through the TRD hits and set XYZ and Pin and Pout
-
-  TTree *hitTree = gAlice->TreeH();
-  gAlice->ResetHits();
-  Int_t   preamp, amp, det, plane = 0, nBytes = 0;
-  Int_t   nTrack = (Int_t) hitTree->GetEntries();
-  Float_t  pos[3], rot[3];
-  Float_t  prepos[3], prerot[3];
-
-  Bool_t  entrance = kFALSE;
-  Bool_t  exit = kFALSE;
-
-
-  Double_t xin[6], xout[6];
-  for(Int_t j = 0; j < 6; j++) {
-    xin[j] = Tracker->GetX(0,j,14);
-    xout[j] = Tracker->GetX(0,j,0);
-  }
-
-  AliTRDhit *hit;
-
-  Int_t current_track = -1;
-
-  printf("\n");
-  while(nTrack--) {
-
-    gAlice->ResetHits();
-    nBytes += hitTree->GetEvent(nTrack);
-    entrance = kTRUE;
-    preamp = 1111;
-    amp = 1111;
-    for(Int_t j = 0; j < 3; j++) { pos[j] = 1111.; rot[j] = 1111.;}
-
-    for(hit = (AliTRDhit *) fTRD->FirstHit(-1); 
-       hit; 
-       hit = (AliTRDhit *) fTRD->NextHit()) {
-
-      preamp = amp;
-      for(Int_t j = 0; j < 3; j++) { prepos[j] = pos[j]; prerot[j] = rot[j];}
-
-      amp = hit->GetCharge();
-      if(amp != 0) continue; 
-      track   = hit->GetTrack(); 
-
-      TParticle *p = gAlice->Particle(track);
-      if (p->Pt() <= low_pt_cut) continue; 
-
-      if(track != current_track) {
-       current_track = track;
-       //      printf("\n\n");
-      }
-
-      det = hit->GetDetector();
-      plane    = fGeo->GetPlane(det); 
-      pos[0] = hit->X();
-      pos[1] = hit->Y();
-      pos[2] = hit->Z();
-      fGeo->Rotate(det,pos,rot);
-
-      if(TMath::Abs(rot[0] - xin[plane]) < 1.2) entrance = kTRUE;
-      else  entrance = kFALSE;
-      if(TMath::Abs(rot[0] - xout[plane]) < 1.2) exit = kTRUE;
-      else  exit = kFALSE;
-      
-      if(entrance && (preamp == 0) && (prerot[0] < 200)) {
-       index = Find(track, TRDmcTracks);
-       if(index > 0) {
-         tpr = 0;
-         tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index);
-         if(tpr) {
-           label = tpr->GetTrackIndex();
-           if(label != track) printf("Got label %d while expected %d !\n",
-                                     label, track);        
-           else {
-             /*
-             printf("\n momentum at plane %d entrance, Pxyz: %f, %f, %f\n",
-                    plane, prepos[0], prepos[1], prepos[2]);
-             printf(" XYZ at plane %d entrance: %f, %f, %f\n",
-                    plane, rot[0], rot[1], rot[2]);
-             printf("expected x range: %f .. %f\n",xin[plane]-1.2,xin[plane]+1.2);
-             */
-             tpr->SetPin(plane, prepos[0], prepos[1], prepos[2]);
-             tpr->SetXYZin(plane, rot[0], rot[1], rot[2]);
-           }
-         }
-       }
-      }
-
-
-      if(exit && (preamp == 0) && (prerot[0] < 200)) {
-       index = Find(track, TRDmcTracks);
-       if(index > 0) {
-         tpr = 0;
-         tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index);
-         if(tpr) {
-           label = tpr->GetTrackIndex();
-           if(label != track) printf("Got label %d while expected %d !\n",
-                                     label, track);        
-           else {
-             /*
-             printf("\n momentum at plane %d exit, Pxyz: %f, %f, %f\n",
-                    plane, prepos[0], prepos[1], prepos[2]);
-             printf(" XYZ at plane %d exit: %f, %f, %f\n",
-                    plane, rot[0], rot[1], rot[2]);
-             printf("expected x range: %f .. %f\n",xout[plane]-1.2,xout[plane]+1.2);
-             */
-             tpr->SetPout(plane, prepos[0], prepos[1], prepos[2]);
-             tpr->SetXYZout(plane, rot[0], rot[1], rot[2]);
-           }
-         }
-       }
-      }
-    }
-  }
-  
-       
-  // Write TRDmcTracks in output file  
-
-  TDirectory *savedir=gDirectory; 
-  Char_t *filename   = "AliTRDmcTracks.root";
-  TFile *out = new TFile(filename,"RECREATE");
-
-  TTree *tracktree = new TTree("MCtracks","TRD MC tracks");
-
-  AliTRDmcTrack *iotrack=0;
-  tracktree->Branch("MCtracks","AliTRDmcTrack",&iotrack,32000,0);
-  
-  Int_t ntracks = TRDmcTracks->GetEntriesFast();
-  
-  for (Int_t i=0; i<ntracks; i++) {
-    AliTRDmcTrack *pt=(AliTRDmcTrack*)TRDmcTracks->UncheckedAt(i);
-    
-    Int_t n = pt->GetNumberOfClusters();
-    if(n > min_no_of_clusters) { 
-      iotrack=pt;
-      tracktree->Fill();
-      printf("Put track with label %d and %d clusters in the output tree \n",
-            pt->GetTrackIndex(),n);
-    }
-  }
-
-  tracktree->Write();
-  out->Close();     
-  savedir->cd(); 
-
-  return;
-}
-
diff --git a/TRD/Macros/AliTRDsaveTrackableSeeds.C b/TRD/Macros/AliTRDsaveTrackableSeeds.C
deleted file mode 100644 (file)
index db5099e..0000000
+++ /dev/null
@@ -1,321 +0,0 @@
-#ifndef __CINT__
-  #include <iostream.h>
-  #include "AliTRDtracker.h"
-
-  #include "AliTRDcluster.h" 
-  #include "AliTRDhit.h" 
-  #include "AliTRDv1.h"
-  #include "AliTRDgeometry.h"    
-  #include "AliTRDparameter.h"    
-  #include "alles.h"  
-  #include "AliTRDmcTrack.h"
-
-  #include "AliTRDtrack.h"
-  #include "AliTrackReference.h"
-  
-  #include "TFile.h"
-  #include "TParticle.h"
-  #include "TStopwatch.h"
-#endif
-
-
-void AliTRDsaveTrackableSeeds() {
-
-  TObjArray mctracks(2000);
-  TObjArray *TRDmcTracks = &mctracks; 
-
-  TFile *geofile =TFile::Open("AliTRDclusters.root");
-  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
-  Int_t nEvent = 0;
-  Tracker->SetEventNumber(nEvent);
-
-  AliTRDgeometry *fGeo   = (AliTRDgeometry*) geofile->Get("TRDgeometry");
-  //AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter");  
-
-  Char_t *alifile = "galice.root";
-
-  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
-  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
-  if (!gafl) {
-    cout << "Open the ALIROOT-file " << alifile << endl;
-    gafl = new TFile(alifile);
-  }
-  else {
-    cout << alifile << " is already open" << endl;
-  }
-
-  // Get AliRun object from file or create it if not on file
-  gAlice = (AliRun*) gafl->Get("gAlice");
-  if (gAlice)
-    cout << "AliRun object found on file" << endl;
-  else
-    gAlice = new AliRun("gAlice","Alice test program");
-
-  AliTRDv1       *fTRD     = (AliTRDv1*) gAlice->GetDetector("TRD");   
-  
-  // Import the Trees for the event nEvent in the file
-  const Int_t nparticles = gAlice->GetEvent(nEvent);
-
-  printf("found %d particles in event %d \n", nparticles, nEvent);
-    
-  // Create TRDmcTracks for each tpc seed
-  Int_t label = -1, charge = 0;
-  Bool_t primary;
-  Float_t mass;
-  Int_t pdg_code;
-
-  TDirectory *savedir=gDirectory;
-
-  TFile *in=TFile::Open("AliTPCBackTracks.root");
-  if (!in->IsOpen()) {
-    cerr<<"can't open file AliTPCBackTracks.root  !\n"; return;
-  }
-
-  char   tname[100];
-  sprintf(tname,"seedsTPCtoTRD_%d",nEvent);
-  TTree *seedTree=(TTree*)in->Get(tname);
-  if (!seedTree) {
-     cerr<<"AliTRDtracker::PropagateBack(): ";
-     cerr<<"can't get a tree with seeds from TPC !\n";
-  }
-
-  AliTPCtrack *seed=new AliTPCtrack;
-  seedTree->SetBranchAddress("tracks",&seed);
-
-  Int_t nSeeds = (Int_t) seedTree->GetEntries();
-  for (Int_t is=0; is<nSeeds; is++) {
-     seedTree->GetEvent(is);
-     Int_t lbl = seed->GetLabel();
-     if(TMath::Abs(lbl) > nparticles) { 
-       printf("extra high seed label %d \n", lbl);
-       continue;
-     }
-     TParticle *p = gAlice->Particle(TMath::Abs(lbl));
-
-     primary = kTRUE;
-     if (p->GetFirstMother()>=0) primary = kFALSE;
-    
-     pdg_code = (Int_t) p->GetPdgCode();
-
-     if ((pdg_code == 10010020) ||
-        (pdg_code == 10010030) || 
-        (pdg_code == 50000050) || 
-        (pdg_code == 50000051) || 
-        (pdg_code == 10020040)) {
-       mass = 0.;
-       charge = 0;
-     }
-     else {
-       TParticlePDG *pdg = p->GetPDG();
-       charge = (Int_t) pdg->Charge(); 
-       mass=pdg->Mass();
-     }
-
-     AliTRDmcTrack *t = new AliTRDmcTrack(TMath::Abs(lbl), lbl, primary, 
-                                         mass, charge, pdg_code);
-     TRDmcTracks->AddLast(t);
-  }
-  delete seed;
-  delete seedTree;
-
-  savedir->cd();                         
-
-  Int_t i, mctIndex[nparticles];
-  for(i=0; i<nparticles; i++) mctIndex[i]=-1;
-
-  Int_t nMCtracks = TRDmcTracks->GetEntriesFast();
-  AliTRDmcTrack *mct = 0;
-  for(i = 0; i < nMCtracks; i++) {
-    mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(i);
-    mctIndex[mct->GetTrackIndex()] = i;
-  }
-
-  // Loop through the TRD rec points and set XYZ and Pin and Pout
-
-  Double_t xin[6], xout[6];
-  Double_t Px, Py, Pz;
-
-  Float_t  pos[3], rot[3];
-
-  for(Int_t j = 0; j < 6; j++) {
-    xin[j] = Tracker->GetX(0,j,0)-3;
-    xout[j] = Tracker->GetX(0,j,0);
-  }
-
-  TTree *trRefTree = gAlice->TreeTR();
-  if (!trRefTree) {
-    cout << "<AliTRDreadTrackRef> No TR tree found" << endl;
-    return;
-  }       
-
-  Int_t nBytes   = 0;
-  Bool_t fEntrance, fExit;
-
-  Int_t nTrack = (Int_t) trRefTree->GetEntries();
-  cout << "<AliTRDreadTrackRef> Found " << nTrack
-       << " primary particles with track refs" << endl;
-
-  for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
-
-    printf("TrackRefs for track %d out of %d \n",iTrack,nTrack);
-    gAlice->ResetTrackReferences();
-    nBytes += trRefTree->GetEvent(iTrack);  
-
-    AliTrackReference *tr = 0;
-    tr = (AliTrackReference*) fTRD->FirstTrackReference(-1);
-    
-    while (tr) {
-
-      Int_t   track = tr->GetTrack();
-      if(mctIndex[track] >= 0) { 
-       mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(mctIndex[track]);
-       
-       pos[0] = tr->X();
-       pos[1] = tr->Y();
-       pos[2] = tr->Z();
-       
-       Double_t phi = TMath::ATan2(pos[1],pos[0]);
-       if(phi < 0) phi = 2 * TMath::Pi() + phi;
-       
-       Int_t sec = ((Int_t) (phi*180/TMath::Pi())) / 20;
-       fGeo->Rotate(fGeo->GetDetector(0,0,17-sec), pos, rot);
-       
-       fEntrance = kFALSE; fExit = kFALSE;
-       for(i = 0; i < 6; i++) {
-         if(TMath::Abs(xin[i] - rot[0]) < 1.4) { fEntrance = kTRUE; break; }
-         if(TMath::Abs(xout[i] - rot[0]) < 1.4) { fExit = kTRUE; break; }
-       }
-      
-       Px = tr->Px();
-       Py = tr->Py();
-       Pz = tr->Pz();
-       
-       if(fEntrance) {
-         mct->SetXYZin(i, rot[0], rot[1], rot[2]);
-         mct->SetPin(i, Px, Py, Pz);
-         /*
-           printf("entr plane %d: rp_x = %f vs %f = xin \n",i,rot[0],xin[i]); 
-           printf("             : Px, Py, Pz = %f, %f, %f \n",Px,Py,Pz); 
-           printf("             : y, z = %f, %f\n",rot[1],rot[2]); 
-         */
-       }
-       if(fExit) {
-         mct->SetXYZout(i, rot[0], rot[1], rot[2]);
-         mct->SetPout(i, Px, Py, Pz);
-         /*
-           printf("exit plane %d: rp_x = %f vs %f = xout \n",i,rot[0],xout[i]); 
-           printf("             : Px, Py, Pz = %f, %f, %f \n",Px,Py,Pz); 
-           printf("             : y, z = %f, %f\n",rot[1],rot[2]);
-         */ 
-       }       
-      }                                   
-      tr = (AliTrackReference *) fTRD->NextTrackReference();
-    }
-  }  
-
-  // Loop through TRD clusters and assign indexes to MC tracks
-  
-  Char_t *clusterfile = "AliTRDclusters.root";
-  TObjArray carray(2000);
-  TObjArray *ClustersArray = &carray;  
-  Tracker->ReadClusters(ClustersArray,clusterfile);
-
-  printf("done with ReadClusters()\n");
-
-  Int_t nClusters = carray.GetEntriesFast();
-
-  Int_t track, det, det0, det1, ltb, plane, ind0, ind1, ncl;
-  Double_t y, y0, y1, q, q0, q1;
-
-  AliTRDcluster *c0 = NULL, *c1 = NULL;
-  printf("nClusters = %d \n", nClusters);
-
-  for (Int_t i = 0; i < nClusters; i++) {
-
-    printf("\r assigning cluster %d out of %d", i, nClusters);
-    AliTRDcluster *cl = (AliTRDcluster *) ClustersArray->UncheckedAt(i);
-
-    for(Int_t j=0; j<3; j++) {
-      track = cl->GetLabel(j);
-      if(track < 0) continue;
-      if(track >= nparticles) { 
-       printf("Track index %d is larger than total number of particles %d\n"
-               ,track,nparticles);
-       continue;
-      }     
-      if(mctIndex[track] < 0) continue;
-      mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(mctIndex[track]);
-
-      label = mct->GetTrackIndex();
-      if(label != track) { 
-       printf("Got label %d while expected %d !\n", label, track);
-       continue;
-      }
-
-      ncl = mct->GetNumberOfClusters();
-      mct->SetNumberOfClusters(ncl+1);
-
-      det=cl->GetDetector();
-      plane = fGeo->GetPlane(det);
-      ltb=cl->GetLocalTimeBin();
-      if((ltb < 0) || (ltb >= kMAX_TB)) continue;
-
-      ind0 = mct->GetClusterIndex(ltb, plane, 0);
-      ind1 = mct->GetClusterIndex(ltb, plane, 1);
-      if(ind0 < 0) {
-       mct->Update(ltb,plane,0,i);            
-      } else {
-       c0 = (AliTRDcluster *) ClustersArray->UncheckedAt(ind0);
-       det0 = c0->GetDetector();
-       y =  cl->GetY();
-       y0 = c0->GetY();
-       q  = TMath::Abs(cl->GetQ());
-       q0 = TMath::Abs(c0->GetQ());
-       if((det == det0) && (TMath::Abs(y-y0) < 1.5)) {
-         if(q > q0) mct->Update(ltb,plane,0,i);
-       } else {
-         if (ind1 < 0) {
-           mct->Update(ltb,plane,1,i);
-         } 
-         else {
-           c1 = (AliTRDcluster *) ClustersArray->UncheckedAt(ind1);
-           det1 = c1->GetDetector();
-           y1 = c1->GetY();
-           q1 = TMath::Abs(c1->GetQ());
-           if((det == det0) && (TMath::Abs(y-y0) < 1.5)) {
-             if(q > q1) mct->Update(ltb,plane,1,i);
-           }     
-         }       
-       }
-      }
-    }
-  }        
-       
-  // Write TRDmcTracks in output file  
-
-  savedir=gDirectory; 
-  Char_t *filename   = "AliTRDtrackableSeeds.root";
-  TFile *out = new TFile(filename,"RECREATE");
-
-  TTree *tracktree = new TTree("MCtracks","TRD MC tracks");
-
-  AliTRDmcTrack *iotrack=0;
-  tracktree->Branch("MCtracks","AliTRDmcTrack",&iotrack,32000,0);
-  
-  Int_t ntracks = TRDmcTracks->GetEntriesFast();
-  
-  for (Int_t i=0; i<ntracks; i++) {
-    AliTRDmcTrack *pt=(AliTRDmcTrack*)TRDmcTracks->UncheckedAt(i);    
-    iotrack=pt;
-    tracktree->Fill();
-    printf("Put track with label %d and %d clusters in the output tree \n",
-          pt->GetTrackIndex(),pt->GetNumberOfClusters());
-  }
-
-  tracktree->Write();
-  out->Close();     
-  savedir->cd(); 
-
-  return;
-}
-
diff --git a/TRD/Macros/AliTRDseedsAnalysis.C b/TRD/Macros/AliTRDseedsAnalysis.C
deleted file mode 100644 (file)
index abb2fd4..0000000
+++ /dev/null
@@ -1,618 +0,0 @@
-#ifndef __CINT__
-  #include <iostream.h>
-  #include "AliTRDtracker.h"
-  #include "AliTRDcluster.h" 
-  #include "AliTRDv1.h"
-  #include "AliTRDgeometry.h"    
-  #include "AliTRDparameter.h"    
-  #include "alles.h"  
-  #include "AliTRDmcTrack.h"
-  #include "AliTRDtrack.h"
-  
-  #include "TFile.h"
-  #include "TParticle.h"
-  #include "TStopwatch.h"
-#endif
-
-
-void AliTRDseedsAnalysis() {
-
-  const Int_t nPrimaries = (Int_t) 3*86030/4;
-  //  const Int_t nPrimaries = 500;
-
-  Bool_t page[6]; for(Int_t i=0; i<6; i++) page[i]=kTRUE;
-
-  const Int_t nPtSteps = 30;
-  const Float_t maxPt = 3.;
-  const Float_t minPt = 0.;
-  
-  const Int_t nEtaSteps = 40;
-  const Float_t maxEta = 1.;
-  const Float_t minEta = -1.;
-  
-  // page[0]
-  TH1F *hNcl=new TH1F("hNcl","Seeds No of Clusters",255,-9.5,500.5); 
-  hNcl->SetFillColor(4);
-  hNcl->SetXTitle("No of clusters"); 
-  hNcl->SetYTitle("counts"); 
-  TH1F *hNtb=new TH1F("hNtb","Seeds No of TB with clusters",160,-9.5,150.5); 
-  hNtb->SetFillColor(4);
-  hNtb->SetXTitle("No of timebins with clusters"); 
-  hNtb->SetYTitle("counts"); 
-  TH1F *hNep=new TH1F("hNep","Seeds end point",160, -9.5, 150.5); 
-  hNep->SetFillColor(4);
-  hNep->SetXTitle("outermost timebin with cluster"); 
-  hNep->SetYTitle("counts"); 
-  TH1F *hNmg=new TH1F("hNmg","Seeds max gap",160, -9.5, 150.5); 
-  hNmg->SetFillColor(4);
-  hNmg->SetXTitle("max gap (tb)"); 
-  hNmg->SetYTitle("counts"); 
-
-  // page[1]
-  Float_t fMin = -5.5, fMax = 134.5;
-  Int_t   iChan = 140; 
-  TH2F *h2ep = new TH2F("h2ep","MC vs RT (end point)",iChan,fMin,fMax,iChan,fMin,fMax); 
-  h2ep->SetMarkerColor(4);
-  h2ep->SetMarkerSize(2);
-  TH2F *h2ntb = new TH2F("h2ntb","MC vs RT (TB with Clusters)",iChan,fMin,fMax,iChan,fMin,fMax); 
-  h2ntb->SetMarkerColor(4);
-  h2ntb->SetMarkerSize(2);
-  TH2F *h2mg = new TH2F("h2mg","MC vs RT (max gap)",iChan/2,fMin,fMax/2,iChan/2,fMin,fMax/2); 
-  h2mg->SetMarkerColor(4);
-  h2mg->SetMarkerSize(2);
-
-  // page[2]
-  TH1F *hPt_all=new TH1F("hPt_all","Seeds Pt",nPtSteps,minPt,maxPt);
-  hPt_all->SetLineColor(4);
-  hPt_all->SetXTitle("Pt (GeV/c)"); 
-  hPt_all->SetYTitle("counts"); 
-
-  TH1F *hPt_short=new TH1F("hPt_short","Short seeds Pt",nPtSteps,minPt,maxPt);
-  hPt_short->SetLineColor(8);
-  TH1F *hPt_long=new TH1F("hPt_long","Long seeds Pt",nPtSteps,minPt,maxPt);
-  hPt_long->SetLineColor(2);
-  
-  TH1F *hrtPt_short=new TH1F("hrtPt_short","RT from short seeds",nPtSteps,minPt,maxPt);
-  hrtPt_short->SetLineColor(8);
-  TH1F *hrtPt_long=new TH1F("hrtPt_long","RT from long seeds",nPtSteps,minPt,maxPt);
-  hrtPt_long->SetLineColor(2);
-
-  TH1F *hEta_all=new TH1F("hEta_all","Seeds Eta",nEtaSteps,minEta,maxEta);
-  hEta_all->SetLineColor(4);
-  hEta_all->SetXTitle("Eta"); 
-  hEta_all->SetYTitle("counts"); 
-  TH1F *hEta_short=new TH1F("hEta_short","Short seeds Eta",nEtaSteps,minEta,maxEta);
-  hEta_short->SetLineColor(8);
-  TH1F *hEta_long=new TH1F("hEta_long","Long seeds Eta",nEtaSteps,minEta,maxEta);
-  hEta_long->SetLineColor(2);
-  
-  TH1F *hrtEta_short=new TH1F("hrtEta_short","RT from short seeds",nEtaSteps,minEta,maxEta);
-  hrtEta_short->SetLineColor(8);
-  TH1F *hrtEta_long=new TH1F("hrtEta_long","RT from long seeds",nEtaSteps,minEta,maxEta);
-  hrtPt_long->SetLineColor(2);
-
-  // page[3]
-  TH1F *hEff_long = new TH1F("hEff_long","Efficiency vs Pt",nPtSteps,minPt,maxPt); 
-  hEff_long->SetLineColor(2); hEff_long->SetLineWidth(2);  
-  hEff_long->SetXTitle("Pt, GeV/c"); 
-  hEff_long->SetYTitle("Efficiency"); 
-  
-  TH1F *hEff_short = new TH1F("hEff_short","Efficiency short",nPtSteps,minPt,maxPt); 
-  hEff_short->SetLineColor(8); hEff_short->SetLineWidth(2);  
-  hEff_short->SetXTitle("Pt, GeV/c"); 
-  hEff_short->SetYTitle("Efficiency"); 
-  
-  TH1F *hEff_long_eta = new TH1F("hEff_long_eta","Efficiency vs Eta",nEtaSteps,minEta,maxEta); 
-  hEff_long_eta->SetLineColor(2); hEff_long_eta->SetLineWidth(2);  
-  hEff_long_eta->SetXTitle("Pt, GeV/c"); 
-  hEff_long_eta->SetYTitle("Efficiency"); 
-  
-  TH1F *hEff_short_eta = new TH1F("hEff_short_eta","Efficiency short",nEtaSteps,minEta,maxEta); 
-  hEff_short_eta->SetLineColor(8); hEff_short_eta->SetLineWidth(2);  
-  hEff_short_eta->SetXTitle("Pt, GeV/c"); 
-  hEff_short_eta->SetYTitle("Efficiency"); 
-  
-  // page[4]
-  TH1F *hY=new TH1F("hY","Y resolution",50,-0.5,0.5);hY->SetLineColor(4);
-  hY->SetLineWidth(2);
-  hY->SetXTitle("(cm)");
-  TH1F *hZ=new TH1F("hZ","Z resolution",80,-4,4);hZ->SetLineColor(4);
-  hZ->SetLineWidth(2);
-  hZ->SetXTitle("(cm)");
-  TH1F *hPhi=new TH1F("hPhi","PHI resolution",100,-20.,20.); 
-  hPhi->SetFillColor(4);
-  hPhi->SetXTitle("(mrad)");
-  TH1F *hLambda=new TH1F("hLambda","Lambda resolution",100,-100,100);
-  hLambda->SetFillColor(17);
-  hLambda->SetXTitle("(mrad)");
-  TH1F *hPt=new TH1F("hPt","Relative Pt resolution",30,-10.,10.);  
-  
-  // page[5]
-  TH1F *hY_n=new TH1F("hY_n","Normalized Y resolution",50,-0.5,0.5); hY_n->SetLineColor(4);
-  hY_n->SetLineWidth(2);
-  hY_n->SetXTitle("  ");
-  TH1F *hZ_n=new TH1F("hZ_n","Normalized Z resolution",50,-0.5,0.5); hZ_n->SetLineColor(2);
-  hZ_n->SetLineWidth(2);
-  hZ_n->SetXTitle("   ");
-  TH1F *hLambda_n=new TH1F("hLambda_n","Normalized Lambda resolution",50,-0.5,0.5);
-  hLambda_n->SetFillColor(17);
-  TH1F *hPt_n=new TH1F("hPt_n","Normalized Pt resolution",50,-0.5,0.5); 
-  
-
-  Int_t nEvent = 0;
-
-  // Load Seeds saved as MC tracks 
-  TObjArray mctarray(2000);
-  TFile *mctf=TFile::Open("AliTRDtrackableSeeds.root");
-  if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDtrackableSeeds.root !\n"; return;}
-  TTree *mctracktree=(TTree*)mctf->Get("MCtracks");
-  TBranch *mctbranch=mctracktree->GetBranch("MCtracks");
-  Int_t const nMCtracks = (Int_t) mctracktree->GetEntries();
-  cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
-  for (Int_t i=0; i<nMCtracks; i++) {
-    AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack();
-    mctbranch->SetAddress(&ioMCtrack);
-    mctracktree->GetEvent(i);
-    mctarray.AddLast(ioMCtrack);
-  }
-  mctf->Close();                 
-
-  TFile *tf=TFile::Open("AliTRDtracks.root");
-
-  if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
-  TObjArray rtarray(2000);
-
-  char   tname[100];
-  sprintf(tname,"TRDb_%d",nEvent);     
-  TTree *tracktree=(TTree*)tf->Get(tname);
-
-  TBranch *tbranch=tracktree->GetBranch("tracks");
-
-  Int_t nRecTracks = (Int_t) tracktree->GetEntries();
-  cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;
-
-  for (Int_t i=0; i<nRecTracks; i++) {
-    AliTRDtrack *iotrack=new AliTRDtrack();
-    tbranch->SetAddress(&iotrack);
-    tracktree->GetEvent(i);
-    rtarray.AddLast(iotrack);
-  }
-  tf->Close();                 
-
-  //  return;
-
-  AliTRDmcTrack *mct;
-  AliTRDtrack *rt;
-  Int_t rtIndex[nMCtracks];
-  for(Int_t j = 0; j < nMCtracks; j++) {
-    mct = (AliTRDmcTrack*) mctarray.UncheckedAt(j);
-    rtIndex[j] = 0;
-    for (Int_t i=0; i<nRecTracks; i++) {
-      rt = (AliTRDtrack*) rtarray.UncheckedAt(i);
-      Int_t label = rt->GetSeedLabel();
-      if(mct->GetTrackIndex() == label) rtIndex[j] = i;
-    }
-  } 
-
-  // Load clusters
-
-  TFile *geofile =TFile::Open("AliTRDclusters.root");   
-  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
-  Tracker->SetEventNumber(nEvent);
-
-  AliTRDgeometry *fGeom   = (AliTRDgeometry*) geofile->Get("TRDgeometry"); 
-  AliTRDparameter *fPar   = (AliTRDparameter*) geofile->Get("TRDparameter"); 
-  
-  Char_t *alifile = "AliTRDclusters.root";
-  TObjArray carray(2000);
-  TObjArray *ClustersArray = &carray;
-  Tracker->ReadClusters(ClustersArray,alifile);   
-
-
-  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
-  alifile = "galice.root";
-  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
-  if (!gafl) {
-    cout << "Open the ALIROOT-file " << alifile << endl;
-    gafl = new TFile(alifile);
-  }
-  else {
-    cout << alifile << " is already open" << endl;
-  }
-  gAlice = (AliRun*) gafl->Get("gAlice");
-  if (gAlice)
-    cout << "AliRun object found on file" << endl;
-  else
-    gAlice = new AliRun("gAlice","Alice test program");
-
-  AliTRDv1       *TRD = (AliTRDv1*) gAlice->GetDetector("TRD");    
-  const Int_t nparticles = gAlice->GetEvent(nEvent);
-  if (nparticles <= 0) return;
-
-  TParticle *p;
-  Bool_t electrons[300000];
-  for(Int_t i = 0; i < 300000; i++) electrons[i] = kFALSE;
-
-  Bool_t mark_electrons = kFALSE;
-  if(mark_electrons) {
-    printf("mark electrons\n");
-
-    for(Int_t i = nPrimaries; i < nparticles; i++) {
-      p = gAlice->Particle(i);
-      if(p->GetMass() > 0.01) continue;
-      if(p->GetMass() < 0.00001) continue;
-      electrons[i] = kTRUE;
-    }
-  }
-
-  AliTRDcluster *cl = 0;
-  Int_t nw, label, index, ti[3];
-  Int_t mct_tbwc, rt_tbwc, mct_max_gap, rt_max_gap, mct_sector, rt_sector; 
-  Int_t mct_max_tb, rt_max_tb, mct_min_tb, rt_min_tb;
-
-  Int_t det, plane, ltb, gtb, gap, max_gap, sector;
-  Double_t Pt, Px, Py, Pz;
-  Double_t mcPt, mcPx, mcPy, mcPz;
-  Double_t x,y,z, mcX;
-  Int_t rtClusters, rtCorrect;
-
-  Int_t nToFind_long, nFound_long, nToFind_short, nFound_short;
-
-  printf("\n");
-
-  Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
-  Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
-  Double_t dx = (Double_t) fPar->GetTimeBinSize();   
-
-  Int_t tbAmp = fPar->GetTimeBefore();
-  Int_t maxAmp = 0;
-  Int_t tbDrift = fPar->GetTimeMax();
-  Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
-
-  tbDrift = TMath::Min(tbDrift,maxDrift);
-  tbAmp = TMath::Min(tbAmp,maxAmp);
-
-  const Int_t nPlanes = fGeom->Nplan();
-  const Int_t tbpp = Tracker->GetTimeBinsPerPlane();
-  const Int_t nTB = tbpp * nPlanes;
-  Int_t mask[nTB][2];
-
-  nToFind_long = 0; nFound_long = 0;
-  nToFind_short = 0; nFound_short = 0;
-
-  for(Int_t i=0; i < nMCtracks; i++) {
-    mct = (AliTRDmcTrack *) mctarray.UncheckedAt(i);
-    label = mct->GetTrackIndex();
-
-    // Waveform of the MC track
-    for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
-    Int_t mct_ncl = mct->GetNumberOfClusters();
-
-    for(ltb = 0; ltb < kMAX_TB; ltb++) {
-      for(plane = 0; plane < nPlanes; plane++) {
-       for(Int_t n = 0; n < 2; n++) {
-         index = mct->GetClusterIndex(ltb,plane,n);
-         if(index < 0) continue; 
-         cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
-
-         for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
-
-         if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
-           printf("mc wrong track label: %d, %d, %d != %d \n",
-                  ti[0], ti[1], ti[2], label);
-         } 
-         det=cl->GetDetector();
-         if(fGeom->GetPlane(det) != plane) 
-           printf("cluster plane = %d != %d expected plane\n", 
-                  fGeom->GetPlane(det), plane);
-         if(cl->GetLocalTimeBin() != ltb) 
-           printf("cluster ltb = %d != %d expected ltb\n", 
-                  cl->GetLocalTimeBin(), ltb);
-         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-         mask[gtb][n] = index;
-       }
-      }
-    }
-    
-    for(plane = 0; plane < nPlanes; plane++) {
-      for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
-       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-       if(mask[gtb][0] > -1) break;
-      }
-      if(ltb >= -tbAmp) break;
-    }  
-    if((plane == nPlanes) && (ltb == -tbAmp-1)) {
-      // printf("warning: for track %d min tb is not found and set to %d!\n",
-      //        label, nTB-1);
-      mct_min_tb = nTB-1;
-      // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
-    }
-    else {
-      mct_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-    }
-
-    for(plane = nPlanes-1 ; plane>=0; plane--) {
-      for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
-       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-       if(mask[gtb][0] > -1) break;
-      }
-      if(ltb < tbDrift) break;
-    }  
-    if((plane == -1) && (ltb == tbDrift)) {
-      //      printf("warning: for track %d max tb is not found and set to 0!\n",label);
-      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
-      mct_max_tb = 0;
-      //      mcX = Tracker->GetX(0,0,tbDrift-1);
-    }
-    else {
-      mct_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
-      //      mcX = Tracker->GetX(0,plane,ltb);
-      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
-      det=cl->GetDetector();
-      sector = fGeom->GetSector(det);  
-      mct_sector = sector;
-    }
-
-    mct_tbwc = 0;
-    mct_max_gap = 0;
-    gap = -1;
-    
-    for(nw = mct_min_tb; nw < mct_max_tb+1; nw++) {
-      gap++;
-      if(mask[nw][0] > -1) {
-       mct_tbwc++;
-       if(gap > mct_max_gap) mct_max_gap = gap;
-       gap = 0;
-      }
-    }
-
-    //  Waveform of the reconstructed track
-    if(rtIndex[i] >= 0) rt = (AliTRDtrack *) rtarray.UncheckedAt(rtIndex[i]);
-
-    for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
-    Int_t rt_ncl = rt->GetNumberOfClusters();
-
-    for(Int_t n = 0; n < rt_ncl; n++) {
-      index = rt->GetClusterIndex(n);
-      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
-      
-      for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
-      
-      if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
-       // printf("rt wrong track label: %d, %d, %d != %d \n", ti[0], ti[1], ti[2], label);
-       continue;
-      } 
-      
-      det=cl->GetDetector();
-      plane = fGeom->GetPlane(det); 
-      ltb = cl->GetLocalTimeBin(); 
-      gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-      mask[gtb][0] = index;
-    }
-
-    for(plane = 0; plane < nPlanes; plane++) {
-      for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
-       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-       if(mask[gtb][0] > -1) break;
-      }
-      if(ltb >= -tbAmp) break;
-    }  
-    if((plane == nPlanes) && (ltb == -tbAmp-1)) {
-      // printf("warning: for track %d min tb is not found and set to %d!\n",
-      //            label, nTB-1);
-      rt_min_tb = nTB-1;
-      // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
-    }
-    else {
-      rt_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-    }
-
-    for(plane = nPlanes-1 ; plane>=0; plane--) {
-      for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
-       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
-       if(mask[gtb][0] > -1) break;
-      }
-      if(ltb < tbDrift) break;
-    }  
-    if((plane == -1) && (ltb == tbDrift)) {
-      // printf("warning: for track %d max tb is not found and set to 0!\n",label);
-      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
-      rt_max_tb = 0;
-      //      mcX = Tracker->GetX(0,0,tbDrift-1);
-    }
-    else {
-      rt_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
-      //      mcX = Tracker->GetX(0,plane,ltb);
-      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
-      det=cl->GetDetector();
-      sector = fGeom->GetSector(det);  
-      rt_sector = sector;
-    }
-
-    rt_tbwc = 0;
-    rt_max_gap = 0;
-    gap = -1;
-    
-    for(nw = rt_min_tb; nw < rt_max_tb+1; nw++) {
-      gap++;
-      if(mask[nw][0] > -1) {
-       rt_tbwc++;
-       if(gap > rt_max_gap) rt_max_gap = gap;
-       gap = 0;
-      }
-    }
-
-    // Fill the histoes
-
-    if(page[0]) {
-      hNcl->Fill((Float_t) mct_ncl);
-      hNtb->Fill((Float_t) mct_tbwc);
-      hNep->Fill((Float_t) mct_max_tb);
-      hNmg->Fill((Float_t) mct_max_gap);
-    }
-    if(page[1]) {
-      h2ep->Fill((Float_t) rt_max_tb, (Float_t) mct_max_tb);
-      h2ntb->Fill((Float_t) rt_tbwc, (Float_t) mct_tbwc);
-      h2mg->Fill((Float_t) rt_max_gap, (Float_t) mct_max_gap);
-    }
-    if(page[2]) {
-      p = gAlice->Particle(label);
-      hPt_all->Fill(p->Pt());
-      hEta_all->Fill(p->Eta());
-      if(mct_max_tb > 60) {
-       nToFind_long++;
-       hPt_long->Fill(p->Pt());
-       hEta_long->Fill(p->Eta());
-       if(((mct_max_tb - rt_max_tb) < 10) && 
-          (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
-         nFound_long++;
-         hrtPt_long->Fill(p->Pt());
-         hrtEta_long->Fill(p->Eta());
-       }         
-      }
-      if((mct_max_tb < 60) && (mct_max_tb > 10)) {
-       nToFind_short++;
-       hPt_short->Fill(p->Pt());
-       hEta_short->Fill(p->Eta());
-       if(((mct_max_tb - rt_max_tb) < 10) && 
-          (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
-         nFound_short++;
-         hrtPt_short->Fill(p->Pt());
-         hrtEta_short->Fill(p->Eta());
-       }         
-      }
-    }
-    if(page[4] && page[5]) {
-      if((mct_tbwc > 50) && (rt_tbwc > 50)) {
-       if(rt->GetSeedLabel() != label) {
-         printf("mct vs rt index mismatch: %d != %d \n",
-                label, rt->GetSeedLabel());
-         return;
-       }
-       Pt = rt->GetPt(); 
-       Double_t cc = TMath::Abs(rt->GetSigmaC2()); 
-       mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1);
-       rt->GetPxPyPz(Px,Py,Pz);      
-
-       printf("\n\ntrack %d \n", label);
-       printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz); 
-       printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz); 
-    
-       mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy);
-       if(mcPt < 0.0001) mcPt = 0.0001;
-    
-       Float_t lamg=TMath::ATan2(mcPz,mcPt);
-       Float_t lam =TMath::ATan2(Pz,Pt);
-       if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
-       else hPt_n->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc));
-    
-       Float_t phig=TMath::ATan2(mcPy,mcPx);
-       Float_t phi =TMath::ATan2(Py,Px);    
-       //      if(!(rt->PropagateTo(x))) continue;
-       //    if((mcPt > Pt_min) && (mcPt < Pt_max)) {
-       hLambda->Fill((lam - lamg)*1000.);
-       hLambda_n->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2()));
-       if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
-       else hPt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.); 
-       hPhi->Fill((phi - phig)*1000.);
-       hY->Fill(rt->GetY() - y);
-       hZ->Fill(rt->GetZ() - z);
-       hY_n->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2()));
-       hZ_n->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2()));
-      }
-    }  
-  }
-
-  // plot pictures
-
-  if(page[0]) {
-    TCanvas *c0=new TCanvas("c0","",0,0,700,850);
-    gStyle->SetOptStat(111110);
-    c0->Divide(2,2);
-    c0->cd(1); gPad->SetLogy(); hNcl->Draw();
-    c0->cd(2); hNtb->Draw();
-    c0->cd(3); hNep->Draw();
-    c0->cd(4); hNmg->Draw();
-    c0->Print("c0.ps","ps"); 
-  }
-  if(page[1]) {
-    TCanvas *c1=new TCanvas("c1","",0,0,700,850);
-    gStyle->SetOptStat(0);
-    c1->Divide(2,2);
-    c1->cd(1); h2ep->Draw();
-    c1->cd(2); h2ntb->Draw();
-    c1->cd(3); h2mg->Draw();
-    //    c1->cd(4); hNmg->Draw();
-    c1->Print("c1.ps","ps"); 
-  }
-  if(page[2]) {
-    TCanvas *c2=new TCanvas("c2","",0,0,700,850);
-    gStyle->SetOptStat(0);
-    c2->Divide(2,2);
-    c2->cd(1); hPt_all->Draw(); hPt_long->Draw("same"); hPt_short->Draw("same"); 
-    c2->cd(2); hEta_all->Draw(); hEta_long->Draw("same"); hEta_short->Draw("same"); 
-    //    c2->cd(3); h2mg->Draw();
-    //    c2->cd(4); hNmg->Draw();
-    c2->Print("c2.ps","ps"); 
-  }
-  if(page[3]) {
-    TCanvas *c3=new TCanvas("c3","",0,0,700,850);
-    gStyle->SetOptStat(0);
-    c3->Divide(1,2);
-    c3->cd(1);
-    hrtPt_long->Sumw2(); hPt_long->Sumw2(); 
-    hEff_long->Divide(hrtPt_long,hPt_long,1,1.,"b");
-    hEff_long->SetMaximum(1.4);
-    hEff_long->SetYTitle("Matching Efficiency");
-    hEff_long->SetXTitle("Pt (GeV/c)");
-    hEff_long->Draw();          
-    hrtPt_short->Sumw2(); hPt_short->Sumw2(); 
-    hEff_short->Divide(hrtPt_short,hPt_short,1,1.,"b");
-    hEff_short->SetMaximum(1.4);
-    hEff_short->SetYTitle("Matching Efficiency");
-    hEff_short->SetXTitle("Pt (GeV/c)");
-    hEff_short->Draw("same");          
-
-    c3->cd(2);
-    hrtEta_long->Sumw2(); hEta_long->Sumw2(); 
-    hEff_long_eta->Divide(hrtEta_long,hEta_long,1,1.,"b");
-    hEff_long_eta->SetMaximum(1.4);
-    hEff_long_eta->SetYTitle("Matching Efficiency");
-    hEff_long_eta->SetXTitle("Eta");
-    hEff_long_eta->Draw();          
-    hrtEta_short->Sumw2(); hEta_short->Sumw2(); 
-    hEff_short_eta->Divide(hrtEta_short,hEta_short,1,1.,"b");
-    hEff_short_eta->SetMaximum(1.4);
-    hEff_short_eta->SetYTitle("Matching Efficiency");
-    hEff_short_eta->SetXTitle("Eta");
-    hEff_short_eta->Draw("same");          
-    c3->Print("c3.ps","ps"); 
-  }
-  if(page[4]) {
-    TCanvas *c4=new TCanvas("c4","",0,0,700,850);
-    gStyle->SetOptStat(111110);
-    c4->Divide(2,3);
-    c4->cd(1); hY->Draw();
-    c4->cd(2); hZ->Draw();
-    c4->cd(3); hPhi->Draw();
-    c4->cd(4); hLambda->Draw();
-    c4->cd(5); hPt->Draw();
-    c4->Print("c4.ps","ps"); 
-  }
-  if(page[5]) {
-    TCanvas *c5=new TCanvas("c5","",0,0,700,850);
-    gStyle->SetOptStat(111110);
-    c5->Divide(2,3);
-    c5->cd(1); hY_n->Draw();
-    c5->cd(2); hZ_n->Draw();
-    c5->cd(3); hLambda_n->Draw();
-    c5->cd(4); hPt_n->Draw();
-    c5->Print("c5.ps","ps"); 
-  }
-  printf("Efficiency (long) = %d/%d = %f\n",nFound_long,nToFind_long,
-        ((Float_t) nFound_long / (Float_t) nToFind_long));
-  printf("Efficiency (shrt) = %d/%d = %f\n",nFound_short,nToFind_short,
-        ((Float_t) nFound_short / (Float_t) nToFind_short));
-}
-
-