[u/mrichter/AliRoot.git] / TPC / TPCDigits2Clusters.C


/*
#include "alles.h"
extern AliTPCParam *gTPCParam;
extern AliTPCClustersArray * gCalcClusters;
extern AliTPCClustersArray * gDifClusters;
AliTPCDigitsArray * GetDigitsArray();
  AliTPCClustersArray * GetExactClustersArray();
AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE);
AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE);
TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters, 
				TClonesArray *diffclusters, Float_t shx, Float_t shy);
*/

void   TPCDigits2Clusters(AliTPCDigitsArray * digarr, AliTPCClustersArray *calcarr,
			  AliTPCClustersArray *exarr=0, AliTPCClustersArray *difarr=0)
{
  //
  //calculate clusters position from digits information
  //

  
  //initialisation of cluster finder
  AliClusterFinder cf;
  cf.SetThreshold(gTPCParam->GetZeroSup());
  cf.SetDetectorParam(gTPCParam);
  cf.SetBFit(kFALSE);
  cf.GetMinuit()->SetPrintLevel(-1);
  cf.GetMinuit()->SetMaxIterations(20);
  
  if (digarr ==0) digarr = GetDigitsArray();
  if ( (digarr == 0) || (digarr->GetTree()==0)) return;  
  if (exarr==0) exarr =GetExactClustersArray();
  if ((exarr!=0)&&(difarr==0)) difarr= GetDifClustersArray();

  if (calcarr==0) calcarr = GetCalcClustersArray(kTRUE);
  //loop over all writen digits segments
  Int_t nsegment = (Int_t)digarr->GetTree()->GetEntries();
  Int_t segment;
  for (segment = 0; segment<nsegment;segment++){
    //load segment with index (TTree internal)  number segment
    AliSimDigits *digrow= (AliSimDigits*)digarr->LoadEntry(segment);
    Int_t sector,row;
    ((AliTPCParam*)digarr->GetParam())->AdjustSectorRow(digrow->GetID(),sector,row);
    
    AliTPCClustersRow * clrow = calcarr->CreateRow(sector,row);
    AliTPCClustersRow * difrow= difarr->CreateRow(sector,row);    
    AliTPCClustersRow * exrow =  exarr->GetRow(sector,row);
    if (exrow==0) exrow = exarr->LoadRow(sector,row);
    TH2F *his = (TH2F*)digrow->GenerHisto();
    digrow->ExpandTrackBuffer();
    if (his==0) return;

    //set current index for cluster finder
    Int_t  index[3]= {0,sector,row};
    cf.GetHisto(his);
    cf.SetDetectorIndex(index);    
    cf.FindPeaks3(clrow->GetArray());

    //get cluster tracks ID
    Int_t ncl = clrow->GetArray()->GetEntriesFast();
    for (Int_t icl = 0 ;icl<ncl; icl++)
      {
	AliDigitCluster *  cl =(AliDigitCluster*)clrow->GetArray()->UncheckedAt(icl);
	Int_t i = (Int_t)cl->fMaxX;
	Int_t j = (Int_t)cl->fMaxY;		
	cl->fTracks[0]=digrow->GetTrackIDFast(i,j,0);
	cl->fTracks[1]=digrow->GetTrackIDFast(i,j,1);
	cl->fTracks[2]=digrow->GetTrackIDFast(i,j,2);	
      }
    Float_t shx = gTPCParam->GetZOffset()/gTPCParam->GetZWidth()+0.5;
    Float_t shy= 0.5;
    CompareClusters(clrow->GetArray(),exrow->GetArray(),difrow->GetArray(),shx,shy);        
    calcarr->StoreRow(sector,row);
    difarr->StoreRow(sector,row);
    digarr->ClearRow(sector,row);
  }    
  //STORING RESULTS
 char treeName[100];
 sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());
 calcarr->GetTree()->Write(treeName);
 if (difarr!=0){
   char treeName[100];
   sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
   difarr->GetTree()->Write(treeName);
 }

}


AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE, const char* name=0) 
{
  //
  //construct AliTPCClusters Array object
  //
  AliTPCClustersArray * arr=0;
  //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {    
  //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
  //  if (arr!=0) arr->Update();
  //}
  if (arr==0) {    
    arr = new AliTPCClustersArray;
    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
    if (file==0) file = TFile::Open("rfio:galice.root","update");   
    arr->SetClusterType("AliDigitCluster");         
    arr->Setup(gTPCParam);
    cout<<"Update status : "<<arr->Update()<<"\n";
    char treeName[100];
    if (name==0)
      sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
    else
      sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());    
    if (newtree!=kTRUE){
      cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
    }
    else {
      arr->MakeTree();
    }
  }
  gCalcClusters = arr;
  return arr;
}


AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE, const char* name=0) 
{
  //
  //construct AliTPCClusters Array object
  //
  AliTPCClustersArray * arr=0;
  //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {    
  //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
  //  if (arr!=0) arr->Update();
  //}
  if (arr==0) {    
    arr = new AliTPCClustersArray;
    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
    if (file==0) file = TFile::Open("rfio:galice.root","update");   
    arr->SetClusterType("AliDifCluster");         
    arr->Setup(gTPCParam);
    cout<<"Update status : "<<arr->Update()<<"\n";
    char treeName[100];
    if (name==0) 
      sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
    else
      sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
    if (newtree!=kTRUE){
      cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
    }
    else {
      arr->MakeTree();
    }
  }
  gDifClusters  = arr; 
  return arr;
}


TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters, 
				TClonesArray *diffclusters, Float_t shx, Float_t shy)
{
  if (calcclusters==0) return 0;
  if (exactclusters==0) return 0;
  if (diffclusters==0) return 0;
  TClonesArray & diff=*diffclusters;
   
   
  Int_t nclusters2 = calcclusters->GetEntriesFast();
  Int_t nclusters1 = exactclusters->GetEntriesFast();
  
  for(Int_t i=0; i<nclusters1; i++){
    AliCluster * cl1 = (AliCluster*)exactclusters->UncheckedAt(i);
    AliDigitCluster *cl2;
    AliDifCluster diffc;
    Int_t index=-1;
    Float_t dx=10000;
    Float_t dy=10000;
    for (Int_t j=0; j<nclusters2;j++){
      cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(j);
      if ( (cl2!=0)&& (cl1!=0)){ 
	Float_t ddx = (cl2->fX-shx)-cl1->fX;
	Float_t ddy = (cl2->fY-shy)-cl1->fY;
	if ((ddx*ddx+ddy*ddy)<(dx*dx+dy*dy)){
	  dx=ddx;
	  dy=ddy;
	  index=j;
	}      
      }
    }
    
    cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(index);
    if (cl2!=0){
      diffc.fDx      =dx;
      diffc.fDy      =dy;
      
      diffc.fAngleX  = cl1->fSigmaX2;
      diffc.fAngleY  = cl1->fSigmaY2; 
      diffc.fTracks[0] = cl2->fTracks[0];
      diffc.fTracks[1] = cl2->fTracks[1];
      diffc.fTracks[2] = cl2->fTracks[2];
      diffc.fGenerTrack= cl1->fTracks[0];

      diffc.fX       = cl2->fX;
      diffc.fY       = cl2->fY;
      diffc.fNx      = cl2->fNx;
      diffc.fNy      = cl2->fNy;
      diffc.fQ       = cl2->fQ;    
      diffc.fOrigQ   = cl1->fQ;
      diffc.fSigmaX2 = cl2->fSigmaX2;
      diffc.fSigmaY2 = cl2->fSigmaY2;
      diffc.fArea    = cl2->fArea;
      diffc.fMax     = cl2->fMax;
    }
    else cout<<"pici[ici/n";
    
    new(diff[i]) AliDifCluster(diffc);
  }
  return diffclusters;
}
Commit	Line	Data
cc80f89e	1
	2	/*
	3	#include "alles.h"
	4	extern AliTPCParam *gTPCParam;
	5	extern AliTPCClustersArray * gCalcClusters;
	6	extern AliTPCClustersArray * gDifClusters;
	7	AliTPCDigitsArray * GetDigitsArray();
	8	AliTPCClustersArray * GetExactClustersArray();
	9	AliTPCClustersArray * GetCalcClustersArray(Bool_t newtree=kFALSE);
	10	AliTPCClustersArray * GetDifClustersArray(Bool_t newtree=kFALSE);
	11	TClonesArray * CompareClusters(TClonesArray calcclusters,TClonesArray exactclusters,
	12	TClonesArray *diffclusters, Float_t shx, Float_t shy);
	13	*/
	14
	15	void TPCDigits2Clusters(AliTPCDigitsArray * digarr, AliTPCClustersArray *calcarr,
	16	AliTPCClustersArray exarr=0, AliTPCClustersArray difarr=0)
	17	{
	18	//
	19	//calculate clusters position from digits information
	20	//
	21
	22
	23	//initialisation of cluster finder
	24	AliClusterFinder cf;
	25	cf.SetThreshold(gTPCParam->GetZeroSup());
	26	cf.SetDetectorParam(gTPCParam);
	27	cf.SetBFit(kFALSE);
	28	cf.GetMinuit()->SetPrintLevel(-1);
	29	cf.GetMinuit()->SetMaxIterations(20);
	30
	31	if (digarr ==0) digarr = GetDigitsArray();
	32	if ( (digarr == 0) \|\| (digarr->GetTree()==0)) return;
	33	if (exarr==0) exarr =GetExactClustersArray();
	34	if ((exarr!=0)&&(difarr==0)) difarr= GetDifClustersArray();
	35
	36	if (calcarr==0) calcarr = GetCalcClustersArray(kTRUE);
	37	//loop over all writen digits segments
	38	Int_t nsegment = (Int_t)digarr->GetTree()->GetEntries();
	39	Int_t segment;
	40	for (segment = 0; segment<nsegment;segment++){
	41	//load segment with index (TTree internal) number segment
	42	AliSimDigits digrow= (AliSimDigits)digarr->LoadEntry(segment);
	43	Int_t sector,row;
	44	((AliTPCParam*)digarr->GetParam())->AdjustSectorRow(digrow->GetID(),sector,row);
	45
	46	AliTPCClustersRow * clrow = calcarr->CreateRow(sector,row);
	47	AliTPCClustersRow * difrow= difarr->CreateRow(sector,row);
	48	AliTPCClustersRow * exrow = exarr->GetRow(sector,row);
	49	if (exrow==0) exrow = exarr->LoadRow(sector,row);
	50	TH2F his = (TH2F)digrow->GenerHisto();
	51	digrow->ExpandTrackBuffer();
	52	if (his==0) return;
	53
	54	//set current index for cluster finder
	55	Int_t index[3]= {0,sector,row};
	56	cf.GetHisto(his);
	57	cf.SetDetectorIndex(index);
	58	cf.FindPeaks3(clrow->GetArray());
	59
	60	//get cluster tracks ID
	61	Int_t ncl = clrow->GetArray()->GetEntriesFast();
	62	for (Int_t icl = 0 ;icl<ncl; icl++)
	63	{
	64	AliDigitCluster * cl =(AliDigitCluster*)clrow->GetArray()->UncheckedAt(icl);
65	Int_t i = (Int_t)cl->fMaxX;
66	Int_t j = (Int_t)cl->fMaxY;
67	cl->fTracks[0]=digrow->GetTrackIDFast(i,j,0);
68	cl->fTracks[1]=digrow->GetTrackIDFast(i,j,1);
69	cl->fTracks[2]=digrow->GetTrackIDFast(i,j,2);
70	}
71	Float_t shx = gTPCParam->GetZOffset()/gTPCParam->GetZWidth()+0.5;
72	Float_t shy= 0.5;
73	CompareClusters(clrow->GetArray(),exrow->GetArray(),difrow->GetArray(),shx,shy);
74	calcarr->StoreRow(sector,row);
75	difarr->StoreRow(sector,row);
76	digarr->ClearRow(sector,row);
77	}
78	//STORING RESULTS
79	char treeName[100];
80	sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());
81	calcarr->GetTree()->Write(treeName);
82	if (difarr!=0){
83	char treeName[100];
84	sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
85	difarr->GetTree()->Write(treeName);
86	}
87
88	}
89
90
91	AliTPCClustersArray * GetCalcClustersArray(Bool_t newtree=kFALSE, const char* name=0)
92	{
93	//
94	//construct AliTPCClusters Array object
95	//
96	AliTPCClustersArray * arr=0;
97	// if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {
98	// arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
99	// if (arr!=0) arr->Update();
100	//}
101	if (arr==0) {
102	arr = new AliTPCClustersArray;
103	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject("galice.root");
106ea0fc	104	if (file==0) file = TFile::Open("rfio:galice.root","update");
cc80f89e	105	arr->SetClusterType("AliDigitCluster");
	106	arr->Setup(gTPCParam);
	107	cout<<"Update status : "<<arr->Update()<<"\n";
	108	char treeName[100];
	109	if (name==0)
	110	sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
	111	else
	112	sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());
	113	if (newtree!=kTRUE){
	114	cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
	115	}
	116	else {
	117	arr->MakeTree();
	118	}
	119	}
	120	gCalcClusters = arr;
	121	return arr;
	122	}
	123
	124
	125	AliTPCClustersArray * GetDifClustersArray(Bool_t newtree=kFALSE, const char* name=0)
	126	{
	127	//
	128	//construct AliTPCClusters Array object
	129	//
	130	AliTPCClustersArray * arr=0;
	131	// if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {
	132	// arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
	133	// if (arr!=0) arr->Update();
	134	//}
	135	if (arr==0) {
	136	arr = new AliTPCClustersArray;
	137	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject("galice.root");
106ea0fc	138	if (file==0) file = TFile::Open("rfio:galice.root","update");
cc80f89e	139	arr->SetClusterType("AliDifCluster");
	140	arr->Setup(gTPCParam);
	141	cout<<"Update status : "<<arr->Update()<<"\n";
	142	char treeName[100];
	143	if (name==0)
	144	sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
	145	else
	146	sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
	147	if (newtree!=kTRUE){
	148	cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
	149	}
	150	else {
	151	arr->MakeTree();
	152	}
	153	}
	154	gDifClusters = arr;
	155	return arr;
	156	}
	157
	158
	159
	160	TClonesArray * CompareClusters(TClonesArray calcclusters,TClonesArray exactclusters,
	161	TClonesArray *diffclusters, Float_t shx, Float_t shy)
	162	{
	163	if (calcclusters==0) return 0;
	164	if (exactclusters==0) return 0;
	165	if (diffclusters==0) return 0;
	166	TClonesArray & diff=*diffclusters;
	167
	168
	169	Int_t nclusters2 = calcclusters->GetEntriesFast();
	170	Int_t nclusters1 = exactclusters->GetEntriesFast();
	171
	172	for(Int_t i=0; i<nclusters1; i++){
	173	AliCluster * cl1 = (AliCluster*)exactclusters->UncheckedAt(i);
	174	AliDigitCluster *cl2;
	175	AliDifCluster diffc;
	176	Int_t index=-1;
	177	Float_t dx=10000;
	178	Float_t dy=10000;
	179	for (Int_t j=0; j<nclusters2;j++){
	180	cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(j);
	181	if ( (cl2!=0)&& (cl1!=0)){
	182	Float_t ddx = (cl2->fX-shx)-cl1->fX;
	183	Float_t ddy = (cl2->fY-shy)-cl1->fY;
	184	if ((ddxddx+ddyddy)<(dxdx+dydy)){
	185	dx=ddx;
	186	dy=ddy;
	187	index=j;
	188	}
	189	}
	190	}
	191
	192	cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(index);
	193	if (cl2!=0){
	194	diffc.fDx =dx;
	195	diffc.fDy =dy;
	196
	197	diffc.fAngleX = cl1->fSigmaX2;
	198	diffc.fAngleY = cl1->fSigmaY2;
	199	diffc.fTracks[0] = cl2->fTracks[0];
	200	diffc.fTracks[1] = cl2->fTracks[1];
	201	diffc.fTracks[2] = cl2->fTracks[2];
	202	diffc.fGenerTrack= cl1->fTracks[0];
203
204	diffc.fX = cl2->fX;
205	diffc.fY = cl2->fY;
206	diffc.fNx = cl2->fNx;
207	diffc.fNy = cl2->fNy;
208	diffc.fQ = cl2->fQ;
209	diffc.fOrigQ = cl1->fQ;
210	diffc.fSigmaX2 = cl2->fSigmaX2;
211	diffc.fSigmaY2 = cl2->fSigmaY2;
212	diffc.fArea = cl2->fArea;
213	diffc.fMax = cl2->fMax;
214	}
215	else cout<<"pici[ici/n";
216
217	new(diff[i]) AliDifCluster(diffc);
218	}
219	return diffclusters;
220	}