ITS/ITStracks.C

   1 #include "AliITSgeom.h"
   2 struct GoodTrack {
   3   Int_t lab;
   4   Int_t code;
   5   Bool_t flag;
   6   Float_t px,py,pz;
   7   Float_t x,y,z;
   8   Float_t pxg,pyg,pzg,ptg;
   9 };
  10
  11
  12 Int_t ITStracks(Int_t evNumber1=0,Int_t evNumber2=0,Int_t nclust=5) {
  13
  14   cerr<<"Select tracks which have nclust rec points in ITS...\n";
  15
  16   GoodTrack gt[15000];
  17   Int_t ngood=0;
  18   ifstream in("good_tracks_tpc");
  19
  20   cerr<<"Reading good tracks...\n";
  21   while (in>>gt[ngood].lab>>gt[ngood].code
  22   >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
  23   >>gt[ngood].x  >>gt[ngood].y >>gt[ngood].z
  24   >>gt[ngood].pxg  >>gt[ngood].pyg >>gt[ngood].pzg
  25   >>gt[ngood].ptg >>gt[ngood].flag) {
  26     ngood++;
  27     cerr<<ngood<<'\r';
  28     if (ngood==15000) {
  29       cerr<<"Too many good tracks !\n";
  30       break;
  31     }
  32   }
  33   if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";
  34
  35
  36   if (gClassTable->GetID("AliRun") < 0) {
  37     gROOT->LoadMacro("loadlibs.C");
  38     loadlibs();
  39   } else {
  40     delete gAlice;
  41     gAlice=0;
  42   }
  43
  44 // Connect the Root Galice file containing Geometry, Kine and Hits
  45   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  46   if (!file) file = new TFile("galice.root","UPDATE");
  47   //if (!file) file = new TFile(filename);
  48
  49 // Get AliRun object from file or create it if not on file
  50   if (!gAlice) {
  51     gAlice = (AliRun*)file->Get("gAlice");
  52     if (gAlice) printf("AliRun object found on file\n");
  53     if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  54   }
  55
  56
  57   AliITS* ITS =(AliITS *)gAlice->GetDetector("ITS");
  58   if (!ITS) return;
  59
  60   //TClonesArray  *recPoints = ITS->RecPoints();
  61   //cout<<" recPoints = "<<recPoints<<"\n";
  62   Int_t numbpoints=0;
  63   AliITSRecPoint *recp=0;
  64   Int_t track=-3;
  65    // Int_t modmin[]={0,80,240,324,500,1282};
  66    // Int_t modmax[]={79,239,323,499,1281,2269};
  67
  68
  69          Int_t modmin[6], modmax[6];
  70          Int_t Nladder[6], Ndetector[6];
  71     AliITSgeom *g1 = ITS->GetITSgeom();
  72
  73     for(Int_t i=0; i<6; i++) {
  74     Nladder[i]=g1->GetNladders(i+1);
  75     Ndetector[i]=g1->GetNdetectors(i+1);
  76          //cout<<Nladder[i]<<"    "<<Ndetector[i]<<"\n";
  77   }
  78      for(Int_t i=0; i<6; i++) {
  79     modmin[i]=g1->GetModuleIndex(i+1,1,1);
  80     modmax[i]=g1->GetModuleIndex(i+1,Nladder[i],Ndetector[i]);
  81          //cout<<modmin[i]<<"    "<<modmax[i]<<"\n";
  82   }
  83
  84
  85
  86   for (int nev=evNumber1; nev<= evNumber2; nev++) {
  87     Int_t nparticles = gAlice->GetEvent(nev);
  88     cout << "nev         " << nev <<endl;
  89     cout << "nparticles  " << nparticles <<endl;
  90     if (nev < evNumber1) continue;
  91     if (nparticles <= 0) return;
  92     TClonesArray  *recPoints = ITS->RecPoints();
  93    TObjArray *particles=gAlice->Particles();
  94    // TClonesArray *particles=gAlice->Particles();
  95    // Int_t np=particles->GetEntriesFast();
  96    Int_t np=gAlice->GetNtrack(); //FCA correction
  97     TMatrix goodITS(np,6);
  98     //Int_t modmin[]={0,80,240,324,500,1282};
  99     //Int_t modmax[]={79,239,323,499,1281,2269};
 100    // for(Int_t j=0; j<np; j++){
 101    //    for(Int_t j1=0; j1<6; j1++) goodITS(j,j1)=0;
 102    // }
 103
 104
 105     Int_t nent=gAlice->TreeR()->GetEntries();
 106     printf("Found %d entries in the TreeR (must be one per module per event!)\n",nent); //
 107
 108     for (Int_t layer=1;layer<=6;layer++) {
 109       for (Int_t mod=modmin[layer-1];mod<=modmax[layer-1];mod++) {
 110              ITS->ResetRecPoints();
 111              //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
 112              gAlice->TreeR()->GetEvent(mod);   //modificato il 2-3-2001
 113              numbpoints = recPoints->GetEntries();
 114              if(!numbpoints) continue;
 115              for (Int_t irec=0;irec<numbpoints;irec++) {
 116                  recp=(AliITSRecPoint*)recPoints->UncheckedAt(irec);
 117                  for (Int_t it=0;it<3;it++) {
 118                      track=recp->GetLabel(it);
 119                      if(track<0) continue;
 120                      if(track > np) {
 121                        cout<<" Error on track number \n";
 122                        getchar();
 123                      }
 124                      goodITS(track,layer-1)=1;
 125                  } //loop over tracks
 126              } //loop over points
 127       } //loop over modules
 128     } //loop over layers
 129
 130     for(Int_t i=0; i<ngood; i++){
 131        if(gt[i].lab>np) {cout<<" Error on TPC good tracks label \n"; getchar();}
 132        Int_t  ilab=gt[i].lab;
 133        Int_t totclust=0;
 134        for(Int_t j=0; j<6; j++) totclust+=goodITS(ilab,j);
 135            if(totclust>=nclust) {
 136              gt[i].flag=1;
 137              printf("label nclusters %d %d\n",gt[i].lab,totclust);  //
 138            }
 139     }
 140
 141     Int_t itsngood=0;
 142     ofstream out("itsgood_tracks");
 143     if (out) {
 144       for (Int_t ngd=0; ngd<ngood; ngd++) {
 145              if(gt[ngd].flag) {
 146                out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
 147                     <<gt[ngd].px <<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
 148                     <<gt[ngd].x  <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<' '
 149                     <<gt[ngd].pxg  <<' '<<gt[ngd].pyg <<' '<<gt[ngd].pzg <<' '
 150                     <<gt[ngd].ptg<<gt[ngd].flag<<endl;
 151                itsngood++;
 152              }
 153       }
 154     } else cerr<<"Can not open file (itsgood_tracks) !\n";
 155     out.close();
 156
 157     cerr<<"Number of good tracks in TPC: "<<ngood<<endl;
 158     cerr<<"Number of good tracks in ITS: "<<itsngood<<endl;
 159
 160   }   // event loop
 161
 162   file->Close();
 163
 164   printf("after file close\n");
 165
 166   return 0;
 167
 168 }
 169