[u/mrichter/AliRoot.git] / ITS / AliITSTracksV1.C

#ifndef __CINT__
#include "AliITSgeom.h"
#include "TParticle.h"
  cerr<<"Reading tp good tracks...\n";
  while (in>>gt[ngood].lab>>gt[ngood].code
  >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
#endif

struct GoodTrack {
  Int_t fEventN;
  Int_t lab;
  Int_t code;
  Bool_t flag;
  Float_t px,py,pz;
  Float_t x,y,z;
  Float_t pxg,pyg,pzg,ptg;
};


Int_t AliITSTracksV1(Int_t evNumber1=0,Int_t evNumber2=0,Int_t nclust=5) {

  cerr<<"Select tracks which have nclust rec points in ITS...\n";

  GoodTrack gt[30000];
  Int_t ngood=0;
  ifstream in("good_tracks_tpc");
   ofstream out("itsgood_tracks");  

  cerr<<"Reading good tracks...\n";
  while (in>>gt[ngood].fEventN>>gt[ngood].lab>>gt[ngood].code
  >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
  >>gt[ngood].x  >>gt[ngood].y >>gt[ngood].z
  >>gt[ngood].pxg  >>gt[ngood].pyg >>gt[ngood].pzg
  >>gt[ngood].ptg >>gt[ngood].flag) {
    ngood++;
    cerr<<ngood<<'\r';
    if (ngood==30000) {
      cerr<<"Too many good tracks !\n";
      break;
    }
  }
  
  if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";

  
   TVector ntrackevtpc(evNumber2-evNumber1 +1);
   Int_t nchange=-1;
   Int_t nmev=-100;
   
   
    Int_t i;
    for( i =0; i<ngood ; i++){ 
    if(gt[i].fEventN != nmev ){
    nmev=gt[i].fEventN;
    nchange++; 
    }
    if(nmev == gt[i].fEventN) ntrackevtpc(nchange)++;
    }
    
 
  if (gClassTable->GetID("AliRun") < 0) {
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
  } else {
    delete gAlice;
    gAlice=0;
  }

// Connect the Root Galice file containing Geometry, Kine and Hits
  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  if (!file) file = new TFile("galice.root","UPDATE");
  //if (!file) file = new TFile(filename);

// Get AliRun object from file or create it if not on file
  if (!gAlice) {
    gAlice = (AliRun*)file->Get("gAlice");
    if (gAlice) printf("AliRun object found on file\n");
    if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  }


  AliITS* ITS =(AliITS *)gAlice->GetDetector("ITS");
  if (!ITS) return;

  //TClonesArray  *recPoints = ITS->RecPoints();
  //cout<<" recPoints = "<<recPoints<<"\n";
  Int_t numbpoints=0;
  AliITSRecPoint *recp=0;
  Int_t track=-3;
   // Int_t modmin[]={0,80,240,324,500,1282};
   // Int_t modmax[]={79,239,323,499,1281,2269};
	 
	
	 Int_t modmin[6], modmax[6];
	 Int_t Nladder[6], Ndetector[6];
    AliITSgeom *g1 = ITS->GetITSgeom(); 
 	 
    for(Int_t i=0; i<6; i++) {
    Nladder[i]=g1->GetNladders(i+1);
    Ndetector[i]=g1->GetNdetectors(i+1);
	 //cout<<Nladder[i]<<"    "<<Ndetector[i]<<"\n";
  }
     for(Int_t i=0; i<6; i++) {
    modmin[i]=g1->GetModuleIndex(i+1,1,1);
    modmax[i]=g1->GetModuleIndex(i+1,Nladder[i],Ndetector[i]);    
	 //cout<<modmin[i]<<"    "<<modmax[i]<<"\n";
  } 
  
  	 
  for (int nev=evNumber1; nev<= evNumber2; nev++) {
    Int_t nparticles = gAlice->GetEvent(nev);
    cout << "nev         " << nev <<endl;
    cout << "nparticles  " << nparticles <<endl;
    if (nev < evNumber1) continue;
    if (nparticles <= 0) return; 
    TClonesArray  *recPoints = ITS->RecPoints();
   TObjArray *particles=gAlice->Particles();     
   // TClonesArray *particles=gAlice->Particles(); 
   // Int_t np=particles->GetEntriesFast();
   Int_t np=gAlice->GetNtrack(); //FCA correction      
    TMatrix goodITS(np,6);
    //Int_t modmin[]={0,80,240,324,500,1282};
    //Int_t modmax[]={79,239,323,499,1281,2269};
   // for(Int_t j=0; j<np; j++){
   //    for(Int_t j1=0; j1<6; j1++) goodITS(j,j1)=0; 
   // }
       
       
    Int_t nent=gAlice->TreeR()->GetEntries();
    printf("Found %d entries in the TreeR (must be one per module per event!)\n",nent); //

    for (Int_t layer=1;layer<=6;layer++) {
      for (Int_t mod=modmin[layer-1];mod<=modmax[layer-1];mod++) {
	     ITS->ResetRecPoints();
	     //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
	     gAlice->TreeR()->GetEvent(mod);   //modificato il 2-3-2001
	     numbpoints = recPoints->GetEntries();
	     if(!numbpoints) continue;
	     for (Int_t irec=0;irec<numbpoints;irec++) {
                 recp=(AliITSRecPoint*)recPoints->UncheckedAt(irec);
		 for (Int_t it=0;it<3;it++) {
		     track=recp->GetLabel(it);
		     if(track<0) continue;
		     if(track > np) { 
		       cout<<" Error on track number \n";
		       getchar();
		     }
		     goodITS(track,layer-1)=1;
	         } //loop over tracks
	     } //loop over points
      } //loop over modules   
    } //loop over layers

    Int_t mingood=0, maxgood=0,jj=0;
    if(nev==evNumber1) mingood=0;
    else
    {for(jj=evNumber1; jj<nev; jj++) mingood += ntrackevtpc(jj);}
     for(jj=evNumber1; jj<=nev; jj++) maxgood+= ntrackevtpc(jj);
    //cout<<" mingood maxgood = "<<mingood<<" "<<maxgood<<"\n";
     for(i=mingood; i<maxgood; i++){
       if(gt[i].lab>np) {cout<<" Error on TPC good tracks label \n"; getchar();}    	   
       Int_t  ilab=gt[i].lab;
       Int_t totclust=0;
       for(Int_t j=0; j<6; j++) totclust+=goodITS(ilab,j);	   
	   if(totclust>=nclust && nev==gt[i].fEventN) {
	     gt[i].flag=1;
	   //  printf("label nclusters %d %d\n",gt[i].lab,totclust);  //
	   }
    }
	   
    Int_t itsngood=0;
 
    if (out) {
     // for (Int_t ngd=0; ngd<ngood; ngd++) {
      for(Int_t ngd=mingood; ngd<maxgood; ngd++){
      TParticle *p = (TParticle*)gAlice->Particle(TMath::Abs(gt[ngd].lab));  //aggiunto 27-9	
		gt[ngd].x=p->Vx(); gt[ngd].y=p->Vy(); 	gt[ngd].z=p->Vz(); 		
	     if(gt[ngd].flag) {
	       out<<gt[ngd].fEventN<<' '<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
		    <<gt[ngd].px <<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
		    <<gt[ngd].x  <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<' '
		    <<gt[ngd].pxg  <<' '<<gt[ngd].pyg <<' '<<gt[ngd].pzg <<' '
		    <<gt[ngd].ptg<<gt[ngd].flag<<endl;
	       itsngood++;
	     }
      }
    } else cerr<<"Can not open file (itsgood_tracks) !\n";


    //cerr<<"Number of good tracks in TPC: "<<ngood<<endl;
    cerr<<"Number of good tracks in TPC: "<<ntrackevtpc(nev)<<endl;
    cerr<<"Number of good tracks in ITS: "<<itsngood<<endl;

  }   // event loop 
  
  out.close();   
  file->Close();   

  printf("after file close\n");
  
  return 0;
  
}
Commit	Line	Data
af8e1c2d	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