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

struct GoodTrack {
  Int_t lab;
  Int_t code;
  Float_t px,py,pz;
  Float_t x,y,z;
  Float_t pxg,pyg,pzg,ptg;
  Bool_t flag;
};
Int_t good_tracks(GoodTrack *gt, Int_t max);


Int_t TPCtracks() {
   Int_t i;
   cerr<<"Doing comparison...\n";

  // Connect the Root Galice file containing Geometry, Kine and Hits


   TFile *cf=TFile::Open("AliTPCclusters.root");
   if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
   AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
   if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }

// Load clusters  
   AliTPCClustersArray *ca=new AliTPCClustersArray;
   ca->Setup(digp);
   ca->SetClusterType("AliTPCcluster");
   ca->ConnectTree("Segment Tree");
   Int_t nentr=Int_t(ca->GetTree()->GetEntries());
   for (Int_t i=0; i<nentr; i++) {
       ca->LoadEntry(i);
   }

// Load tracks
   TFile *tf=TFile::Open("AliTPCtracks.root");
   if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
   TObjArray tarray(2000);
   TTree *tracktree=(TTree*)tf->Get("TreeT");
   TBranch *tbranch=tracktree->GetBranch("tracks");
   nentr=(Int_t)tracktree->GetEntries();
   for (i=0; i<nentr; i++) {
       AliTPCtrack *iotrack=new AliTPCtrack;
       tbranch->SetAddress(&iotrack);
       tracktree->GetEvent(i);
       iotrack->CookLabel(ca);
       tarray.AddLast(iotrack);
   }   

   tf->Close();
   delete ca;
   cf->Close();
   //printf("after cf close\n");
   
   cerr<<"Number of found tracks "<<nentr<<endl;
   tarray.Sort();

   TFile *pfile = new TFile("tpctracks.root","RECREATE"); 

   TTree tracktree1("TreeT","Tree with TPC tracks");
   AliTPCtrack *iotrack=0;   
   tracktree1.Branch("tracks","AliTPCtrack",&iotrack,32000,0);

   for (i=0; i<nentr; i++) {
       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
       tracktree1.Fill();
   }   

   tracktree1.Write();
   pfile->Close();
   tarray.Clear();


/////////////////////////////////////////////////////////////////////////
   GoodTrack gt[15000];
   Int_t ngood=0;

   cerr<<"Marking good tracks (this will take a while)...\n";
   ngood=good_tracks(gt,15000);
   ofstream out("good_tracks");
   if (out) {
         for (Int_t ngd=0; ngd<ngood; ngd++)            
	    out<<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;
   } else cerr<<"Can not open file (good_tracks) !\n";
   out.close();

   cerr<<"Number of good tracks : "<<ngood<<endl;

   printf("before return in tpctracks\n");
   return 0;
}

//---------------------------------

Int_t good_tracks(GoodTrack *gt, Int_t max) {
   Int_t nt=0;

   //TFile *file=TFile::Open("rfio:galice.root");  // via server
   TFile *file=TFile::Open("galice.root");
   if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}

   if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
     cerr<<"gAlice have not been found on galice.root !\n";
     exit(5);
   }

   gAlice->GetEvent(0);   

   AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
   Int_t ver = TPC->IsVersion(); 
   cerr<<"TPC version "<<ver<<" has been found !\n";

   AliTPCParam *digp=(AliTPCParam*)file->Get("75x40_100x60");
   if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
   TPC->SetParam(digp);

   Int_t nrow_up=digp->GetNRowUp();
   Int_t nrows=digp->GetNRowLow()+nrow_up;
   Int_t zero=digp->GetZeroSup();
   Int_t gap=Int_t(0.125*nrows);
   Int_t good_number=Int_t(0.4*nrows);

   TObjArray *particles=gAlice->Particles(); 
   //Int_t np=particles->GetEntriesFast();
   Int_t np=gAlice->GetNtrack(); //FCA correction   
   Int_t *good=new Int_t[np];
   for (Int_t ii=0; ii<np; ii++) good[ii]=0;
   
   //MI change to be possible compile macro
   //definition out of the switch statement
   Int_t sectors_by_rows=0;
   TTree *TD=0;
   AliSimDigits da, *digits=&da;
   Int_t *count=0;

   switch (ver) {
   case 1:
     {
      TFile *cf=TFile::Open("AliTPCclusters.root");
      if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
      AliTPCClustersArray *ca=new AliTPCClustersArray;
      ca->Setup(digp);
      ca->SetClusterType("AliTPCcluster");
      ca->ConnectTree("Segment Tree");
      Int_t nrows=Int_t(ca->GetTree()->GetEntries());
      for (Int_t n=0; n<nrows; n++) {
          AliSegmentID *s=ca->LoadEntry(n);
          Int_t sec,row;
          digp->AdjustSectorRow(s->GetID(),sec,row);
          AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
          Int_t ncl=clrow.GetArray()->GetEntriesFast();
          while (ncl--) {
              AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
              Int_t lab=c->GetLabel(0);
              if (lab<0) continue; //noise cluster
              lab=TMath::Abs(lab);
              if (sec>=digp->GetNInnerSector())
              if (row==nrow_up-1    ) good[lab]|=0x1000;
              if (sec>=digp->GetNInnerSector())
              if (row==nrow_up-1-gap) good[lab]|=0x800;
              good[lab]++;
          }
          ca->ClearRow(sec,row);
      }
      cf->Close();
     }
      break;
   case 2:
      TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
      TD->GetBranch("Segment")->SetAddress(&digits);
      count = new Int_t[np];
      Int_t i;
      for (i=0; i<np; i++) count[i]=0;
      Int_t sectors_by_rows=(Int_t)TD->GetEntries();
      for (i=0; i<sectors_by_rows; i++) {
          if (!TD->GetEvent(i)) continue;
          Int_t sec,row;
          digp->AdjustSectorRow(digits->GetID(),sec,row);
          cerr<<sec<<' '<<row<<"                                     \r";
          digits->First();
          while (digits->Next()) {
              Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
              Short_t dig = digits->GetDigit(it,ip);
              Int_t idx0=digits->GetTrackID(it,ip,0); 
              Int_t idx1=digits->GetTrackID(it,ip,1);
              Int_t idx2=digits->GetTrackID(it,ip,2);
              if (idx0>=0 && dig>=zero) count[idx0]+=1;
              if (idx1>=0 && dig>=zero) count[idx1]+=1;
              if (idx2>=0 && dig>=zero) count[idx2]+=1;
          }
          for (Int_t j=0; j<np; j++) {
              if (count[j]>1) {
                 if (sec>=digp->GetNInnerSector())
                 if (row==nrow_up-1    ) good[j]|=0x1000;
                 if (sec>=digp->GetNInnerSector())
                 if (row==nrow_up-1-gap) good[j]|=0x800;
                 good[j]++;
              }
              count[j]=0;
          }
      }
      delete[] count;
      break;
   default:
      cerr<<"Invalid TPC version !\n";
      file->Close();
      exit(7);
   }

   TTree *TH=gAlice->TreeH();
   //TClonesArray *hits=TPC->Hits();
   Int_t npart=TH->GetEntries();

   while (npart--) {
      AliTPChit *hit0=0;
      TPC->ResetHits();
      TH->GetEvent(npart);
      AliTPChit *hit = (AliTPChit*) TPC->FirstHit(-1);
      
      while(hit) {
         if(hit->fQ==0.) break;
	 hit = (AliTPChit*) TPC->NextHit();
      }
      if(hit) {
         hit0 = new AliTPChit(*hit); // Make copy of hit
	 hit=hit0;
      }
      else continue;
      AliTPChit *hit1=(AliTPChit*) TPC->NextHit();
      if(hit1==0) continue;
      
      if (hit1->fQ != 0.) continue;
     // Int_t i=hit->fTrack;
	   Int_t i=hit->Track();  
      //TParticle *p = (TParticle*)particles->UncheckedAt(i);
	  TParticle *p = (TParticle*)gAlice->Particle(i);  //mod. secondonew
	                                                      //AliTPCcomaprison 
      if (p->GetFirstMother()>=0) continue;  //secondary particle
      if (good[i] < 0x1000+0x800+2+good_number) continue;
      if (p->Pt()<0.100) continue;
      if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;

      gt[nt].lab=i;
      gt[nt].code=p->GetPdgCode();
//**** px py pz - in global coordinate system, x y z - in local !
     // gt[nt].px=hit->fX; gt[nt].py=hit->fY; gt[nt].pz=hit->fZ;  //modificato tenendo conto di AliTPCComparison
      gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();	  
      Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
      //gt[nt].x = hit1->fX*cs + hit1->fY*sn;
     // gt[nt].y =-hit1->fX*sn + hit1->fY*cs;  //modificato tenedo conto di AliTPCComaprison
      //gt[nt].z = hit1->fZ;
      gt[nt].x = hit1->X()*cs + hit1->Y()*sn;
      gt[nt].y =-hit1->X()*sn + hit1->Y()*cs;
      gt[nt].z = hit1->Z();		
      gt[nt].pxg = p->Px();
      gt[nt].pyg = p->Py();
      gt[nt].pzg = p->Pz();				
      gt[nt].ptg = p->Pt();		
      gt[nt].flag = 0;		
      nt++;
        
	if(hit0) delete hit0;
      cerr<<i<<"                \r";
      if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
   }
   delete[] good;

   printf("before delete gAlice\n");

   delete gAlice; gAlice=0;

   printf("after delete gAlice\n");
   file->Close();
   printf("after file close\n");
   return nt;
}

//--------------------------------------
Commit	Line	Data
d1c69928	1	struct GoodTrack {
	2	Int_t lab;
	3	Int_t code;
	4	Float_t px,py,pz;
	5	Float_t x,y,z;
	6	Float_t pxg,pyg,pzg,ptg;
	7	Bool_t flag;
	8	};
	9	Int_t good_tracks(GoodTrack *gt, Int_t max);
	10
	11
	12	Int_t TPCtracks() {
cc6e75dd	13	Int_t i;
d1c69928	14	cerr<<"Doing comparison...\n";
	15
	16	// Connect the Root Galice file containing Geometry, Kine and Hits
	17
	18
d1c69928	19	TFile *cf=TFile::Open("AliTPCclusters.root");
	20	if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
	21	AliTPCParam digp= (AliTPCParam)cf->Get("75x40_100x60");
	22	if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }
	23
	24	// Load clusters
	25	AliTPCClustersArray *ca=new AliTPCClustersArray;
	26	ca->Setup(digp);
	27	ca->SetClusterType("AliTPCcluster");
	28	ca->ConnectTree("Segment Tree");
	29	Int_t nentr=Int_t(ca->GetTree()->GetEntries());
	30	for (Int_t i=0; i<nentr; i++) {
cc6e75dd	31	ca->LoadEntry(i);
d1c69928	32	}
	33
	34	// Load tracks
	35	TFile *tf=TFile::Open("AliTPCtracks.root");
	36	if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
	37	TObjArray tarray(2000);
	38	TTree tracktree=(TTree)tf->Get("TreeT");
	39	TBranch *tbranch=tracktree->GetBranch("tracks");
cc6e75dd	40	nentr=(Int_t)tracktree->GetEntries();
d1c69928	41	for (i=0; i<nentr; i++) {
	42	AliTPCtrack *iotrack=new AliTPCtrack;
	43	tbranch->SetAddress(&iotrack);
	44	tracktree->GetEvent(i);
	45	iotrack->CookLabel(ca);
	46	tarray.AddLast(iotrack);
	47	}
	48
36d58ccf	49	tf->Close();
	50	delete ca;
	51	cf->Close();
	52	//printf("after cf close\n");
	53
d1c69928	54	cerr<<"Number of found tracks "<<nentr<<endl;
	55	tarray.Sort();
	56
	57	TFile *pfile = new TFile("tpctracks.root","RECREATE");
	58
	59	TTree tracktree1("TreeT","Tree with TPC tracks");
36d58ccf	60	AliTPCtrack *iotrack=0;
d1c69928	61	tracktree1.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
	62
	63	for (i=0; i<nentr; i++) {
	64	iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
	65	tracktree1.Fill();
	66	}
	67
	68	tracktree1.Write();
d1c69928	69	pfile->Close();
d1c69928	70	tarray.Clear();
d1c69928	71
d1c69928	72
d1c69928	73	/////////////////////////////////////////////////////////////////////////
2b7ebd74	74	GoodTrack gt[15000];
d1c69928	75	Int_t ngood=0;
	76
	77	cerr<<"Marking good tracks (this will take a while)...\n";
2b7ebd74	78	ngood=good_tracks(gt,15000);
d1c69928	79	ofstream out("good_tracks");
	80	if (out) {
	81	for (Int_t ngd=0; ngd<ngood; ngd++)
	82	out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
	83	<<gt[ngd].px <<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
	84	<<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<' '
	85	<<gt[ngd].pxg <<' '<<gt[ngd].pyg <<' '<<gt[ngd].pzg <<' '
	86	<<gt[ngd].ptg <<gt[ngd].flag<<endl;
	87	} else cerr<<"Can not open file (good_tracks) !\n";
	88	out.close();
	89
	90	cerr<<"Number of good tracks : "<<ngood<<endl;
	91
	92	printf("before return in tpctracks\n");
d1c69928	93	return 0;
	94	}
	95
	96	//---------------------------------
	97
	98	Int_t good_tracks(GoodTrack *gt, Int_t max) {
	99	Int_t nt=0;
	100
	101	//TFile *file=TFile::Open("rfio:galice.root"); // via server
	102	TFile *file=TFile::Open("galice.root");
	103	if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
	104
	105	if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
	106	cerr<<"gAlice have not been found on galice.root !\n";
	107	exit(5);
	108	}
	109
	110	gAlice->GetEvent(0);
	111
	112	AliTPC TPC=(AliTPC)gAlice->GetDetector("TPC");
	113	Int_t ver = TPC->IsVersion();
	114	cerr<<"TPC version "<<ver<<" has been found !\n";
	115
	116	AliTPCParam digp=(AliTPCParam)file->Get("75x40_100x60");
	117	if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
	118	TPC->SetParam(digp);
	119
	120	Int_t nrow_up=digp->GetNRowUp();
	121	Int_t nrows=digp->GetNRowLow()+nrow_up;
	122	Int_t zero=digp->GetZeroSup();
	123	Int_t gap=Int_t(0.125*nrows);
	124	Int_t good_number=Int_t(0.4*nrows);
	125
2b7ebd74	126	TObjArray *particles=gAlice->Particles();
	127	//Int_t np=particles->GetEntriesFast();
	128	Int_t np=gAlice->GetNtrack(); //FCA correction
d1c69928	129	Int_t *good=new Int_t[np];
d1c69928	130	for (Int_t ii=0; ii<np; ii++) good[ii]=0;
cc6e75dd	131
	132	//MI change to be possible compile macro
	133	//definition out of the switch statement
	134	Int_t sectors_by_rows=0;
	135	TTree *TD=0;
	136	AliSimDigits da, *digits=&da;
	137	Int_t *count=0;
d1c69928	138
	139	switch (ver) {
	140	case 1:
	141	{
	142	TFile *cf=TFile::Open("AliTPCclusters.root");
	143	if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
	144	AliTPCClustersArray *ca=new AliTPCClustersArray;
	145	ca->Setup(digp);
	146	ca->SetClusterType("AliTPCcluster");
	147	ca->ConnectTree("Segment Tree");
	148	Int_t nrows=Int_t(ca->GetTree()->GetEntries());
	149	for (Int_t n=0; n<nrows; n++) {
	150	AliSegmentID *s=ca->LoadEntry(n);
	151	Int_t sec,row;
	152	digp->AdjustSectorRow(s->GetID(),sec,row);
	153	AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
	154	Int_t ncl=clrow.GetArray()->GetEntriesFast();
	155	while (ncl--) {
	156	AliTPCcluster c=(AliTPCcluster)clrow[ncl];
	157	Int_t lab=c->GetLabel(0);
	158	if (lab<0) continue; //noise cluster
	159	lab=TMath::Abs(lab);
	160	if (sec>=digp->GetNInnerSector())
	161	if (row==nrow_up-1 ) good[lab]\|=0x1000;
	162	if (sec>=digp->GetNInnerSector())
	163	if (row==nrow_up-1-gap) good[lab]\|=0x800;
	164	good[lab]++;
	165	}
	166	ca->ClearRow(sec,row);
	167	}
	168	cf->Close();
	169	}
	170	break;
	171	case 2:
cc6e75dd	172	TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
d1c69928	173	TD->GetBranch("Segment")->SetAddress(&digits);
cc6e75dd	174	count = new Int_t[np];
d1c69928	175	Int_t i;
	176	for (i=0; i<np; i++) count[i]=0;
	177	Int_t sectors_by_rows=(Int_t)TD->GetEntries();
	178	for (i=0; i<sectors_by_rows; i++) {
	179	if (!TD->GetEvent(i)) continue;
	180	Int_t sec,row;
	181	digp->AdjustSectorRow(digits->GetID(),sec,row);
	182	cerr<<sec<<' '<<row<<" \r";
	183	digits->First();
	184	while (digits->Next()) {
	185	Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
	186	Short_t dig = digits->GetDigit(it,ip);
	187	Int_t idx0=digits->GetTrackID(it,ip,0);
	188	Int_t idx1=digits->GetTrackID(it,ip,1);
	189	Int_t idx2=digits->GetTrackID(it,ip,2);
	190	if (idx0>=0 && dig>=zero) count[idx0]+=1;
	191	if (idx1>=0 && dig>=zero) count[idx1]+=1;
	192	if (idx2>=0 && dig>=zero) count[idx2]+=1;
	193	}
	194	for (Int_t j=0; j<np; j++) {
	195	if (count[j]>1) {
	196	if (sec>=digp->GetNInnerSector())
	197	if (row==nrow_up-1 ) good[j]\|=0x1000;
	198	if (sec>=digp->GetNInnerSector())
	199	if (row==nrow_up-1-gap) good[j]\|=0x800;
	200	good[j]++;
	201	}
	202	count[j]=0;
	203	}
	204	}
	205	delete[] count;
	206	break;
	207	default:
	208	cerr<<"Invalid TPC version !\n";
	209	file->Close();
	210	exit(7);
	211	}
	212
	213	TTree *TH=gAlice->TreeH();
cc6e75dd	214	//TClonesArray *hits=TPC->Hits();
d1c69928	215	Int_t npart=TH->GetEntries();
	216
	217	while (npart--) {
cc6e75dd	218	AliTPChit *hit0=0;
d1c69928	219	TPC->ResetHits();
d1c69928	220	TH->GetEvent(npart);
cc6e75dd	221	AliTPChit hit = (AliTPChit) TPC->FirstHit(-1);
	222
	223	while(hit) {
	224	if(hit->fQ==0.) break;
	225	hit = (AliTPChit*) TPC->NextHit();
	226	}
	227	if(hit) {
	228	hit0 = new AliTPChit(*hit); // Make copy of hit
	229	hit=hit0;
d1c69928	230	}
cc6e75dd	231	else continue;
	232	AliTPChit hit1=(AliTPChit) TPC->NextHit();
	233	if(hit1==0) continue;
	234
d1c69928	235	if (hit1->fQ != 0.) continue;
d1c69928	236	// Int_t i=hit->fTrack;
2b7ebd74	237	Int_t i=hit->Track();
	238	//TParticle p = (TParticle)particles->UncheckedAt(i);
	239	TParticle p = (TParticle)gAlice->Particle(i); //mod. secondonew
	240	//AliTPCcomaprison
d1c69928	241	if (p->GetFirstMother()>=0) continue; //secondary particle
	242	if (good[i] < 0x1000+0x800+2+good_number) continue;
	243	if (p->Pt()<0.100) continue;
	244	if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
	245
	246	gt[nt].lab=i;
	247	gt[nt].code=p->GetPdgCode();
	248	//**** px py pz - in global coordinate system, x y z - in local !
	249	// gt[nt].px=hit->fX; gt[nt].py=hit->fY; gt[nt].pz=hit->fZ; //modificato tenendo conto di AliTPCComparison
	250	gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();
	251	Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
	252	//gt[nt].x = hit1->fXcs + hit1->fYsn;
	253	// gt[nt].y =-hit1->fXsn + hit1->fYcs; //modificato tenedo conto di AliTPCComaprison
	254	//gt[nt].z = hit1->fZ;
	255	gt[nt].x = hit1->X()cs + hit1->Y()sn;
	256	gt[nt].y =-hit1->X()sn + hit1->Y()cs;
	257	gt[nt].z = hit1->Z();
	258	gt[nt].pxg = p->Px();
	259	gt[nt].pyg = p->Py();
	260	gt[nt].pzg = p->Pz();
	261	gt[nt].ptg = p->Pt();
	262	gt[nt].flag = 0;
	263	nt++;
36d58ccf	264
36d58ccf	265	if(hit0) delete hit0;
d1c69928	266	cerr<<i<<" \r";
	267	if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
	268	}
	269	delete[] good;
	270
	271	printf("before delete gAlice\n");
	272
	273	delete gAlice; gAlice=0;
	274
	275	printf("after delete gAlice\n");
	276	file->Close();
	277	printf("after file close\n");
	278	return nt;
	279	}
	280
	281	//--------------------------------------
	282