TRD/Macros/AliTRDtoTOFanalysis.C

   1 #ifndef __CINT__
   2   #include <iostream.h>
   3
   4   #include "AliTOF.h"
   5   #include "AliTOFhit.h"
   6   #include "AliTPCtrack.h"
   7   #include "AliTRDtrack.h"
   8   #include "AliTRDgeometry.h"
   9
  10   #include "alles.h"
  11   #include "TFile.h"
  12   #include "TParticle.h"
  13   #include "TStopwatch.h"
  14 #endif
  15
  16 void AliTRDtoTOFanalysis()
  17 {
  18
  19   Int_t nEvent = 0;
  20   const Int_t maxIndex = 80000;  // max number of primaries to be analysed
  21   Int_t rtIndex[maxIndex];
  22   for(Int_t i = 0; i < maxIndex; i++) rtIndex[i] = -1;
  23
  24   //****** tracking:  Load reconstructed tracks
  25
  26   TFile *tf=TFile::Open("AliTRDtracks.root");
  27
  28   if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
  29   TObjArray tarray(2000);
  30
  31   char   tname[100];
  32   sprintf(tname,"seedsTRDtoTOF2_%d",nEvent);
  33   //  for other radial positions use different trees:
  34   //    sprintf(tname,"seedsTRDtoTOF1_%d",nEvent);
  35   //    sprintf(tname,"seedsTRDtoPHOS_%d",nEvent);
  36   //    sprintf(tname,"seedsTRDtoRHIC_%d",nEvent);
  37
  38   TTree *tracktree=(TTree*)tf->Get(tname);
  39
  40   TBranch *tbranch=tracktree->GetBranch("tracks");
  41
  42   Int_t nRecTracks = (Int_t) tracktree->GetEntries();
  43   cout<<"Found "<<nRecTracks<<" entries in the track tree "<<tname<<endl;
  44
  45   for (Int_t i=0; i<nRecTracks; i++) {
  46     AliTPCtrack *iotrack=new AliTPCtrack();
  47     tbranch->SetAddress(&iotrack);
  48     tracktree->GetEvent(i);
  49     tarray.AddLast(iotrack);
  50     Int_t trackLabel = iotrack->GetLabel();
  51
  52     //    printf("rt with %d clusters and label %d \n",
  53     //     iotrack->GetNumberOfClusters(), trackLabel);
  54
  55     if(trackLabel < 0) continue;
  56     if(trackLabel >= maxIndex) continue;
  57     rtIndex[trackLabel] = i;
  58   }
  59   tf->Close();
  60
  61   //********
  62
  63   TH1F *hdx = new TH1F("dx","X(hit) - X(track)",40,-10,10);
  64   TH1F *hdy = new TH1F("dy","Y(hit) - Y(track)",100,-20,20);
  65   TH1F *hdz = new TH1F("dz","Z(hit) - Z(track)",100,-20,20);
  66   TH2F *hdt = new TH2F("dt","T(hit) vs S(track)",100,0,500,500,200,700);
  67
  68
  69   // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
  70   Char_t *alifile = "galice.root";
  71
  72   TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
  73   if (!gafl) {
  74     cout << "Open the ALIROOT-file " << alifile << endl;
  75     gafl = new TFile(alifile);
  76   }
  77   else {
  78     cout << alifile << " is already open" << endl;
  79   }
  80
  81   // Get AliRun object from file or create it if not on file
  82   gAlice = (AliRun*) gafl->Get("gAlice");
  83   if (gAlice)
  84     cout << "AliRun object found on file" << endl;
  85   else
  86     gAlice = new AliRun("gAlice","Alice test program");
  87
  88
  89   // Import the Trees for the event nEvent in the file
  90   const Int_t nparticles = gAlice->GetEvent(nEvent);
  91   if (nparticles <= 0) return;
  92
  93    Float_t x,y,z,mass,e;
  94    Double_t gX, gY, gZ;   // rec. track coordinates in the global system
  95    Double_t Px, Py, Pz;   // rec. track momentum components in the global system
  96    Float_t X, Y, Z;     // local tracking coordinates for reconstructed track
  97    Int_t nbytes = 0;
  98    Int_t j,hit,ipart;
  99    Int_t nhits;
 100    Float_t tof;
 101    TParticle *particle;
 102    AliTPCtrack *rt;
 103
 104 // Get pointers to Alice detectors and Hits containers
 105    AliDetector *TOF  = gAlice->GetDetector("TOF");
 106    Int_t ntracks    = (Int_t) gAlice->TreeH()->GetEntries();
 107
 108 // Start loop on tracks in the hits containers
 109
 110    for (Int_t track=0; track < ntracks; track++) {
 111
 112      if(TOF) {
 113        for(AliTOFhit* tofHit = (AliTOFhit*)TOF->FirstHit(track);
 114            tofHit;
 115            tofHit=(AliTOFhit*)TOF->NextHit()) {
 116
 117          ipart    = tofHit->GetTrack();
 118          if(ipart >= maxIndex) continue;
 119          if(rtIndex[ipart] < 0) continue;
 120
 121          x = tofHit->X();
 122          y = tofHit->Y();
 123          z = tofHit->Z();
 124
 125          //******* tracking:  extract track coordinates, momentum, etc.
 126
 127          rt = (AliTPCtrack*) tarray.UncheckedAt(rtIndex[ipart]);
 128          Double_t tr_length = rt->GetLength();  // meaningfull only with ITS
 129                                                 // part fully implemented
 130
 131          AliTRDtrack *trd_rt = new AliTRDtrack(*rt, rt->GetAlpha());
 132          trd_rt->GetGlobalXYZ(gX,gY,gZ);  // returns running global coordinates
 133          trd_rt->GetPxPyPz(Px,Py,Pz);  // returns momentum in global system
 134
 135          /*
 136          printf("\n hit position - rec. track position:\n");
 137          printf(" X: %f - %f = %f \n", x, gX, x - gX);
 138          printf(" Y: %f - %f = %f \n", y, gY, y - gY);
 139          printf(" Z: %f - %f = %f \n", z, gZ, z - gZ);
 140          printf("\n");
 141          */
 142
 143          delete trd_rt;
 144
 145          // conversion from hit coordinates from global coordinate system
 146          // to "local" tracking system
 147
 148          Float_t phi =TMath::ATan2(y,x);
 149          if (phi > 2.*TMath::Pi()) phi -= 2.*TMath::Pi();
 150          if (phi < 0.            ) phi += 2.*TMath::Pi();
 151          Int_t sec = Int_t(phi/AliTRDgeometry::GetAlpha()) % 18;
 152          Float_t alpha = (sec + 0.5) * AliTRDgeometry::GetAlpha();
 153
 154          Float_t tmp=x*TMath::Cos(alpha) + y*TMath::Sin(alpha);
 155          y= - x*TMath::Sin(alpha) + y*TMath::Cos(alpha);
 156          x=tmp;
 157
 158          //         particle = (TParticle*)gAlice->Particle(ipart);
 159          //        if (particle->GetFirstMother() < 0) continue;
 160
 161          X = (Float_t) rt->GetX();   // radial position in the tracking system
 162          Y = (Float_t) rt->GetY();   // r*phi position within the sector
 163          Z = (Float_t) rt->GetZ();   // Z position
 164
 165
 166          if(TMath::Abs(X-x-1) > 2) continue;
 167
 168          hdx->Fill(X-x);
 169          hdy->Fill(Y-y);
 170          hdz->Fill(Z-z);
 171
 172          //**********
 173
 174        }
 175      }
 176    }
 177
 178   TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940);
 179   c1->SetFillColor(10);
 180   c1->Divide(2,2);
 181   c1->cd(1); hdx->Draw();
 182   c1->cd(2); hdy->Draw();
 183   c1->cd(3); hdz->Draw();
 184   c1->cd(4); hdt->Draw();
 185 }