ITS/AliITSTracksV1.C

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