Additional include path (Linux icc)
[u/mrichter/AliRoot.git] / 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