[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[7000];
   Int_t ngood=0;

   cerr<<"Marking good tracks (this will take a while)...\n";
   ngood=good_tracks(gt,7000);
   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);

   TClonesArray *particles=gAlice->Particles(); 
   Int_t np=particles->GetEntriesFast();
   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();  //modificata in accordo a nuovo AliTPCComparison
      TParticle *p = (TParticle*)particles->UncheckedAt(i);
      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
	73	/////////////////////////////////////////////////////////////////////////
	74	GoodTrack gt[7000];
	75	Int_t ngood=0;
	76
	77	cerr<<"Marking good tracks (this will take a while)...\n";
	78	ngood=good_tracks(gt,7000);
	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
	126	TClonesArray *particles=gAlice->Particles();
	127	Int_t np=particles->GetEntriesFast();
	128	Int_t *good=new Int_t[np];
	129	for (Int_t ii=0; ii<np; ii++) good[ii]=0;
cc6e75dd	130
	131	//MI change to be possible compile macro
	132	//definition out of the switch statement
	133	Int_t sectors_by_rows=0;
	134	TTree *TD=0;
	135	AliSimDigits da, *digits=&da;
	136	Int_t *count=0;
d1c69928	137
	138	switch (ver) {
	139	case 1:
	140	{
	141	TFile *cf=TFile::Open("AliTPCclusters.root");
	142	if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
	143	AliTPCClustersArray *ca=new AliTPCClustersArray;
	144	ca->Setup(digp);
	145	ca->SetClusterType("AliTPCcluster");
	146	ca->ConnectTree("Segment Tree");
	147	Int_t nrows=Int_t(ca->GetTree()->GetEntries());
	148	for (Int_t n=0; n<nrows; n++) {
	149	AliSegmentID *s=ca->LoadEntry(n);
	150	Int_t sec,row;
	151	digp->AdjustSectorRow(s->GetID(),sec,row);
	152	AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
	153	Int_t ncl=clrow.GetArray()->GetEntriesFast();
	154	while (ncl--) {
	155	AliTPCcluster c=(AliTPCcluster)clrow[ncl];
	156	Int_t lab=c->GetLabel(0);
	157	if (lab<0) continue; //noise cluster
	158	lab=TMath::Abs(lab);
	159	if (sec>=digp->GetNInnerSector())
	160	if (row==nrow_up-1 ) good[lab]\|=0x1000;
	161	if (sec>=digp->GetNInnerSector())
	162	if (row==nrow_up-1-gap) good[lab]\|=0x800;
	163	good[lab]++;
	164	}
	165	ca->ClearRow(sec,row);
	166	}
	167	cf->Close();
	168	}
	169	break;
	170	case 2:
cc6e75dd	171	TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
d1c69928	172	TD->GetBranch("Segment")->SetAddress(&digits);
cc6e75dd	173	count = new Int_t[np];
d1c69928	174	Int_t i;
	175	for (i=0; i<np; i++) count[i]=0;
	176	Int_t sectors_by_rows=(Int_t)TD->GetEntries();
	177	for (i=0; i<sectors_by_rows; i++) {
	178	if (!TD->GetEvent(i)) continue;
	179	Int_t sec,row;
	180	digp->AdjustSectorRow(digits->GetID(),sec,row);
	181	cerr<<sec<<' '<<row<<" \r";
	182	digits->First();
	183	while (digits->Next()) {
	184	Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
	185	Short_t dig = digits->GetDigit(it,ip);
	186	Int_t idx0=digits->GetTrackID(it,ip,0);
	187	Int_t idx1=digits->GetTrackID(it,ip,1);
	188	Int_t idx2=digits->GetTrackID(it,ip,2);
	189	if (idx0>=0 && dig>=zero) count[idx0]+=1;
	190	if (idx1>=0 && dig>=zero) count[idx1]+=1;
	191	if (idx2>=0 && dig>=zero) count[idx2]+=1;
	192	}
	193	for (Int_t j=0; j<np; j++) {
	194	if (count[j]>1) {
	195	if (sec>=digp->GetNInnerSector())
	196	if (row==nrow_up-1 ) good[j]\|=0x1000;
	197	if (sec>=digp->GetNInnerSector())
	198	if (row==nrow_up-1-gap) good[j]\|=0x800;
	199	good[j]++;
	200	}
	201	count[j]=0;
	202	}
	203	}
	204	delete[] count;
	205	break;
	206	default:
	207	cerr<<"Invalid TPC version !\n";
	208	file->Close();
	209	exit(7);
	210	}
	211
	212	TTree *TH=gAlice->TreeH();
cc6e75dd	213	//TClonesArray *hits=TPC->Hits();
d1c69928	214	Int_t npart=TH->GetEntries();
	215
	216	while (npart--) {
cc6e75dd	217	AliTPChit *hit0=0;
d1c69928	218	TPC->ResetHits();
d1c69928	219	TH->GetEvent(npart);
cc6e75dd	220	AliTPChit hit = (AliTPChit) TPC->FirstHit(-1);
	221
	222	while(hit) {
	223	if(hit->fQ==0.) break;
	224	hit = (AliTPChit*) TPC->NextHit();
	225	}
	226	if(hit) {
	227	hit0 = new AliTPChit(*hit); // Make copy of hit
	228	hit=hit0;
d1c69928	229	}
cc6e75dd	230	else continue;
	231	AliTPChit hit1=(AliTPChit) TPC->NextHit();
	232	if(hit1==0) continue;
	233
d1c69928	234	if (hit1->fQ != 0.) continue;
	235	// Int_t i=hit->fTrack;
	236	Int_t i=hit->Track(); //modificata in accordo a nuovo AliTPCComparison
	237	TParticle p = (TParticle)particles->UncheckedAt(i);
	238	if (p->GetFirstMother()>=0) continue; //secondary particle
	239	if (good[i] < 0x1000+0x800+2+good_number) continue;
	240	if (p->Pt()<0.100) continue;
	241	if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
	242
	243	gt[nt].lab=i;
	244	gt[nt].code=p->GetPdgCode();
	245	//**** px py pz - in global coordinate system, x y z - in local !
	246	// gt[nt].px=hit->fX; gt[nt].py=hit->fY; gt[nt].pz=hit->fZ; //modificato tenendo conto di AliTPCComparison
	247	gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();
	248	Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
	249	//gt[nt].x = hit1->fXcs + hit1->fYsn;
	250	// gt[nt].y =-hit1->fXsn + hit1->fYcs; //modificato tenedo conto di AliTPCComaprison
	251	//gt[nt].z = hit1->fZ;
	252	gt[nt].x = hit1->X()cs + hit1->Y()sn;
	253	gt[nt].y =-hit1->X()sn + hit1->Y()cs;
	254	gt[nt].z = hit1->Z();
	255	gt[nt].pxg = p->Px();
	256	gt[nt].pyg = p->Py();
	257	gt[nt].pzg = p->Pz();
	258	gt[nt].ptg = p->Pt();
	259	gt[nt].flag = 0;
	260	nt++;
36d58ccf	261
36d58ccf	262	if(hit0) delete hit0;
d1c69928	263	cerr<<i<<" \r";
	264	if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
	265	}
	266	delete[] good;
	267
	268	printf("before delete gAlice\n");
	269
	270	delete gAlice; gAlice=0;
	271
	272	printf("after delete gAlice\n");
	273	file->Close();
	274	printf("after file close\n");
	275	return nt;
	276	}
	277
	278	//--------------------------------------
	279