[u/mrichter/AliRoot.git] / STEER / STEER / AliMillePedeRecord.cxx

#include "AliMillePedeRecord.h"
#include <TMath.h>
#include "AliLog.h"

/**********************************************************************************************/
/* AliMillePedeRecords: class to store the data of single track processing                    */
/* Format: for each measured point the data is stored consequtively                           */
/* INDEX                                                      VALUE                           */
/* -1                                                         residual                        */
/* Local_param_id                                             dResidual/dLocal_param          */
/* ...                                                        ...                             */
/* -2                                                         weight of the measurement       */
/* Global_param_od                                            dResidual/dGlobal_param         */
/* ...                                                        ...                             */
/*                                                                                            */
/* The records for all processed tracks are stored in the temporary tree in orgder to be      */
/* reused for multiple iterations of MillePede                                                */
/*                                                                                            */ 
/* Author: ruben.shahoyan@cern.ch                                                             */
/*                                                                                            */ 
/**********************************************************************************************/

ClassImp(AliMillePedeRecord)

//_____________________________________________________________________________________________
AliMillePedeRecord::AliMillePedeRecord() : 
fSize(0),fNGroups(0),fRunID(0),fGroupID(0),fIndex(0),fValue(0),fWeight(1) {SetUniqueID(0);}

//_____________________________________________________________________________________________
AliMillePedeRecord::AliMillePedeRecord(const AliMillePedeRecord& src) : 
  TObject(src),fSize(src.fSize),fNGroups(src.fNGroups),fRunID(src.fRunID),fGroupID(0),fIndex(0),fValue(0),fWeight(src.fWeight)
{
  // copy ct-r
  fIndex = new Int_t[GetDtBufferSize()];
  memcpy(fIndex,src.fIndex,fSize*sizeof(Int_t));
  fValue = new Double_t[GetDtBufferSize()];
  memcpy(fValue,src.fValue,fSize*sizeof(Double_t));
  fGroupID = new UShort_t[GetGrBufferSize()];
  memcpy(fGroupID,src.fGroupID,GetGrBufferSize()*sizeof(UShort_t));
}

//_____________________________________________________________________________________________
AliMillePedeRecord& AliMillePedeRecord::operator=(const AliMillePedeRecord& rhs)
{ 
  // assignment op-r
  if (this!=&rhs) {
    Reset();
    for (int i=0;i<rhs.GetSize();i++) {
      Double_t val;
      Int_t    ind;
      rhs.GetIndexValue(i,ind,val);
      AddIndexValue(ind,val);
    }
    fWeight = rhs.fWeight;
    fRunID = rhs.fRunID;
    for (int i=0;i<rhs.GetNGroups();i++) MarkGroup(rhs.GetGroupID(i));
  }
  return *this;
}

//_____________________________________________________________________________________________
AliMillePedeRecord::~AliMillePedeRecord() {delete[] fIndex; delete[] fValue; delete[] fGroupID;}

//_____________________________________________________________________________________________
void AliMillePedeRecord::Reset()
{
  // reset all
  fSize = 0;
  for (int i=fNGroups;i--;) fGroupID[i] = 0;
  fNGroups = 0;
  fRunID = 0;
  fWeight = 1.;
}

//_____________________________________________________________________________________________
void AliMillePedeRecord::Print(const Option_t *) const
{
  // print itself
  if (!fSize) {AliInfo("No data"); return;}
  int cnt=0,point=0;
  //  
  if (fNGroups) printf("Groups: ");
  for (int i=0;i<fNGroups;i++) printf("%4d |",GetGroupID(i)); 
  printf("Run: %9d Weight: %+.2e\n",fRunID,fWeight);
  while(cnt<fSize) {
    //
    Double_t resid = fValue[cnt++];
    Double_t *derLoc = GetValue()+cnt;
    int    *indLoc = GetIndex()+cnt;
    int     nLoc = 0;
    while(!IsWeight(cnt)) {nLoc++;cnt++;}
    Double_t weight = GetValue(cnt++);
    Double_t *derGlo = GetValue()+cnt;
    int    *indGlo = GetIndex()+cnt;
    int     nGlo = 0;
    while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;} 
    //
    printf("\n*** Point#%2d | Residual = %+.4e | Weight = %+.4e\n",point++,resid,weight);
    printf("Locals : "); 
    for (int i=0;i<nLoc;i++) printf("[%5d] %+.4e|",indLoc[i],derLoc[i]); printf("\n");
    printf("Globals: "); 
    for (int i=0;i<nGlo;i++) printf("[%5d] %+.4e|",indGlo[i],derGlo[i]); printf("\n");
    //
  }
  //
}

//_____________________________________________________________________________________________
Double_t AliMillePedeRecord::GetGloResWProd(Int_t indx) const
{
  // get sum of derivative over global variable indx * res. at point * weight
  if (!fSize) {AliInfo("No data"); return 0;}
  int cnt=0;
  double prodsum = 0.0;
  //  
  while(cnt<fSize) {
    //
    Double_t resid = fValue[cnt++];
    while(!IsWeight(cnt)) cnt++;
    Double_t weight = GetValue(cnt++);
    Double_t *derGlo = GetValue()+cnt;
    int    *indGlo = GetIndex()+cnt;
    int     nGlo = 0;
    while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;} 
    for (int i=nGlo;i--;) if (indGlo[i]==indx) prodsum += resid*weight*derGlo[i];
    //
  }
  return prodsum;
}

//_____________________________________________________________________________________________
Double_t AliMillePedeRecord::GetGlobalDeriv(Int_t pnt, Int_t indx) const
{
  // get derivative over global variable indx at point pnt
  if (!fSize) {AliError("No data"); return 0;}
  int cnt=0,point=0;
  //  
  while(cnt<fSize) {
    //
    cnt++;
    while(!IsWeight(cnt)) cnt++;
    cnt++;
    Double_t *derGlo = GetValue()+cnt;
    int    *indGlo = GetIndex()+cnt;
    int     nGlo = 0;
    while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;} 
    //
    if (pnt != point++) continue;
    for (int i=nGlo;i--;) if (indGlo[i]==indx) return derGlo[i];
    break;
  }
  return 0;
  //
}

//_____________________________________________________________________________________________
Double_t AliMillePedeRecord::GetLocalDeriv(Int_t pnt, Int_t indx) const
{
  // get derivative over local variable indx at point pnt
  if (!fSize) {AliError("No data"); return 0;}
  int cnt=0,point=0;
  //  
  while(cnt<fSize) {
    //
    cnt++;
    Double_t *derLoc = GetValue()+cnt;
    int    *indLoc = GetIndex()+cnt;
    int     nLoc = 0;
    while(!IsWeight(cnt)) {nLoc++;cnt++;}
    cnt++;
    while(!IsResidual(cnt) && cnt<fSize) cnt++;
    if (pnt != point++) continue;
    for (int i=nLoc;i--;) if (indLoc[i]==indx) return derLoc[i];
    break;
  }
  return 0;
  //
}

//_____________________________________________________________________________________________
Double_t AliMillePedeRecord::GetResidual(Int_t pnt) const
{
  // get residual at point pnt
  if (!fSize) {AliError("No data"); return 0;}
  int cnt=0,point=0;
  //  
  while(cnt<fSize) {
    //
    Double_t resid = fValue[cnt++];
    while(!IsWeight(cnt)) cnt++;
    cnt++;
    while(!IsResidual(cnt) && cnt<fSize) cnt++;
    if (pnt != point++) continue;
    return resid;
  }
  return 0;
  //
}

//_____________________________________________________________________________________________
Double_t AliMillePedeRecord::GetWeight(Int_t pnt) const
{
  // get weight of point pnt
  if (!fSize) {AliError("No data"); return 0;}
  int cnt=0,point=0;
  //  
  while(cnt<fSize) {
    //
    cnt++;
    while(!IsWeight(cnt)) cnt++;
    if (point==pnt) return GetValue(cnt);;
    cnt++;
    while(!IsResidual(cnt) && cnt<fSize) cnt++;
    point++;
  }
  return -1;
  //
}

//_____________________________________________________________________________________________
void AliMillePedeRecord::ExpandDtBuffer(Int_t bfsize)
{
  // add extra space for derivatives data
  bfsize = TMath::Max(bfsize,GetDtBufferSize());
  Int_t *tmpI = new Int_t[bfsize];
  memcpy(tmpI,fIndex, fSize*sizeof(Int_t));
  delete [] fIndex;
  fIndex = tmpI;
  //
  Double_t *tmpD = new Double_t[bfsize];
  memcpy(tmpD,fValue, fSize*sizeof(Double_t));
  delete [] fValue;
  fValue = tmpD;
  //
  SetDtBufferSize(bfsize);
}

//_____________________________________________________________________________________________
void AliMillePedeRecord::ExpandGrBuffer(Int_t bfsize)
{
  // add extra space for groupID data 
  bfsize = TMath::Max(bfsize,GetGrBufferSize());
  UShort_t *tmpI = new UShort_t[bfsize];
  memcpy(tmpI,fGroupID, fNGroups*sizeof(UShort_t));
  delete [] fGroupID;
  fGroupID = tmpI;
  for (int i=fNGroups;i<bfsize;i++) fGroupID[i] = 0;
  //
  SetGrBufferSize(bfsize);
}

//_____________________________________________________________________________________________
void AliMillePedeRecord::MarkGroup(Int_t id)
{
  // mark the presence of the detector group
  id++; // groupID is stored as realID+1
  if (fNGroups>0 && fGroupID[fNGroups-1]==id) return; // already there
  if (fNGroups>=GetGrBufferSize()) ExpandGrBuffer(2*(fNGroups+1));
  fGroupID[fNGroups++] = id;  
}
Commit	Line	Data
de34b538	1	#include "AliMillePedeRecord.h"
	2	#include <TMath.h>
	3	#include "AliLog.h"
	4
	5	/**********************************************************************************************/
	6	/* AliMillePedeRecords: class to store the data of single track processing */
	7	/* Format: for each measured point the data is stored consequtively */
	8	/* INDEX VALUE */
	9	/* -1 residual */
	10	/* Local_param_id dResidual/dLocal_param */
	11	/* ... ... */
	12	/* -2 weight of the measurement */
	13	/* Global_param_od dResidual/dGlobal_param */
	14	/* ... ... */
	15	/* */
	16	/* The records for all processed tracks are stored in the temporary tree in orgder to be */
	17	/* reused for multiple iterations of MillePede */
	18	/* */
	19	/* Author: ruben.shahoyan@cern.ch */
	20	/* */
	21	/**********************************************************************************************/
	22
	23	ClassImp(AliMillePedeRecord)
	24
	25	//_____________________________________________________________________________________________
	26	AliMillePedeRecord::AliMillePedeRecord() :
e30a812f	27	fSize(0),fNGroups(0),fRunID(0),fGroupID(0),fIndex(0),fValue(0),fWeight(1) {SetUniqueID(0);}
de34b538	28
	29	//_____________________________________________________________________________________________
	30	AliMillePedeRecord::AliMillePedeRecord(const AliMillePedeRecord& src) :
3b94b44f	31	TObject(src),fSize(src.fSize),fNGroups(src.fNGroups),fRunID(src.fRunID),fGroupID(0),fIndex(0),fValue(0),fWeight(src.fWeight)
de34b538	32	{
339fbe23	33	// copy ct-r
de34b538	34	fIndex = new Int_t[GetDtBufferSize()];
	35	memcpy(fIndex,src.fIndex,fSize*sizeof(Int_t));
	36	fValue = new Double_t[GetDtBufferSize()];
	37	memcpy(fValue,src.fValue,fSize*sizeof(Double_t));
a8c5b94c	38	fGroupID = new UShort_t[GetGrBufferSize()];
a8c5b94c	39	memcpy(fGroupID,src.fGroupID,GetGrBufferSize()*sizeof(UShort_t));
de34b538	40	}
	41
	42	//_____________________________________________________________________________________________
	43	AliMillePedeRecord& AliMillePedeRecord::operator=(const AliMillePedeRecord& rhs)
339fbe23	44	{
339fbe23	45	// assignment op-r
de34b538	46	if (this!=&rhs) {
	47	Reset();
	48	for (int i=0;i<rhs.GetSize();i++) {
	49	Double_t val;
	50	Int_t ind;
	51	rhs.GetIndexValue(i,ind,val);
	52	AddIndexValue(ind,val);
	53	}
a8c5b94c	54	fWeight = rhs.fWeight;
e30a812f	55	fRunID = rhs.fRunID;
de34b538	56	for (int i=0;i<rhs.GetNGroups();i++) MarkGroup(rhs.GetGroupID(i));
	57	}
	58	return *this;
	59	}
	60
	61	//_____________________________________________________________________________________________
	62	AliMillePedeRecord::~AliMillePedeRecord() {delete[] fIndex; delete[] fValue; delete[] fGroupID;}
	63
	64	//_____________________________________________________________________________________________
	65	void AliMillePedeRecord::Reset()
	66	{
339fbe23	67	// reset all
de34b538	68	fSize = 0;
a8c5b94c	69	for (int i=fNGroups;i--;) fGroupID[i] = 0;
de34b538	70	fNGroups = 0;
e30a812f	71	fRunID = 0;
a8c5b94c	72	fWeight = 1.;
de34b538	73	}
	74
	75	//_____________________________________________________________________________________________
	76	void AliMillePedeRecord::Print(const Option_t *) const
	77	{
339fbe23	78	// print itself
de34b538	79	if (!fSize) {AliInfo("No data"); return;}
	80	int cnt=0,point=0;
	81	//
	82	if (fNGroups) printf("Groups: ");
a8c5b94c	83	for (int i=0;i<fNGroups;i++) printf("%4d \|",GetGroupID(i));
e30a812f	84	printf("Run: %9d Weight: %+.2e\n",fRunID,fWeight);
de34b538	85	while(cnt<fSize) {
	86	//
	87	Double_t resid = fValue[cnt++];
	88	Double_t *derLoc = GetValue()+cnt;
	89	int *indLoc = GetIndex()+cnt;
	90	int nLoc = 0;
	91	while(!IsWeight(cnt)) {nLoc++;cnt++;}
	92	Double_t weight = GetValue(cnt++);
	93	Double_t *derGlo = GetValue()+cnt;
	94	int *indGlo = GetIndex()+cnt;
	95	int nGlo = 0;
	96	while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;}
	97	//
	98	printf("\n*** Point#%2d \| Residual = %+.4e \| Weight = %+.4e\n",point++,resid,weight);
	99	printf("Locals : ");
	100	for (int i=0;i<nLoc;i++) printf("[%5d] %+.4e\|",indLoc[i],derLoc[i]); printf("\n");
	101	printf("Globals: ");
	102	for (int i=0;i<nGlo;i++) printf("[%5d] %+.4e\|",indGlo[i],derGlo[i]); printf("\n");
	103	//
	104	}
	105	//
	106	}
	107
9b86a63d	108	//_____________________________________________________________________________________________
	109	Double_t AliMillePedeRecord::GetGloResWProd(Int_t indx) const
	110	{
	111	// get sum of derivative over global variable indx * res. at point * weight
	112	if (!fSize) {AliInfo("No data"); return 0;}
	113	int cnt=0;
	114	double prodsum = 0.0;
	115	//
	116	while(cnt<fSize) {
	117	//
	118	Double_t resid = fValue[cnt++];
	119	while(!IsWeight(cnt)) cnt++;
	120	Double_t weight = GetValue(cnt++);
	121	Double_t *derGlo = GetValue()+cnt;
	122	int *indGlo = GetIndex()+cnt;
	123	int nGlo = 0;
	124	while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;}
	125	for (int i=nGlo;i--;) if (indGlo[i]==indx) prodsum += residweightderGlo[i];
	126	//
	127	}
	128	return prodsum;
	129	}
	130
	131	//_____________________________________________________________________________________________
	132	Double_t AliMillePedeRecord::GetGlobalDeriv(Int_t pnt, Int_t indx) const
	133	{
	134	// get derivative over global variable indx at point pnt
	135	if (!fSize) {AliError("No data"); return 0;}
	136	int cnt=0,point=0;
	137	//
	138	while(cnt<fSize) {
	139	//
	140	cnt++;
	141	while(!IsWeight(cnt)) cnt++;
	142	cnt++;
	143	Double_t *derGlo = GetValue()+cnt;
	144	int *indGlo = GetIndex()+cnt;
	145	int nGlo = 0;
	146	while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;}
	147	//
	148	if (pnt != point++) continue;
	149	for (int i=nGlo;i--;) if (indGlo[i]==indx) return derGlo[i];
	150	break;
	151	}
	152	return 0;
	153	//
	154	}
	155
	156	//_____________________________________________________________________________________________
	157	Double_t AliMillePedeRecord::GetLocalDeriv(Int_t pnt, Int_t indx) const
	158	{
	159	// get derivative over local variable indx at point pnt
	160	if (!fSize) {AliError("No data"); return 0;}
	161	int cnt=0,point=0;
	162	//
	163	while(cnt<fSize) {
	164	//
	165	cnt++;
	166	Double_t *derLoc = GetValue()+cnt;
	167	int *indLoc = GetIndex()+cnt;
	168	int nLoc = 0;
	169	while(!IsWeight(cnt)) {nLoc++;cnt++;}
	170	cnt++;
	171	while(!IsResidual(cnt) && cnt<fSize) cnt++;
172	if (pnt != point++) continue;
173	for (int i=nLoc;i--;) if (indLoc[i]==indx) return derLoc[i];
174	break;
175	}
176	return 0;
177	//
178	}
179
180	//_____________________________________________________________________________________________
181	Double_t AliMillePedeRecord::GetResidual(Int_t pnt) const
182	{
183	// get residual at point pnt
184	if (!fSize) {AliError("No data"); return 0;}
185	int cnt=0,point=0;
186	//
187	while(cnt<fSize) {
188	//
189	Double_t resid = fValue[cnt++];
190	while(!IsWeight(cnt)) cnt++;
191	cnt++;
192	while(!IsResidual(cnt) && cnt<fSize) cnt++;
193	if (pnt != point++) continue;
194	return resid;
195	}
196	return 0;
197	//
198	}
199
ca7f226b	200	//_____________________________________________________________________________________________
	201	Double_t AliMillePedeRecord::GetWeight(Int_t pnt) const
	202	{
	203	// get weight of point pnt
	204	if (!fSize) {AliError("No data"); return 0;}
	205	int cnt=0,point=0;
	206	//
	207	while(cnt<fSize) {
	208	//
f2e76498	209	cnt++;
ca7f226b	210	while(!IsWeight(cnt)) cnt++;
	211	if (point==pnt) return GetValue(cnt);;
	212	cnt++;
	213	while(!IsResidual(cnt) && cnt<fSize) cnt++;
	214	point++;
	215	}
	216	return -1;
	217	//
	218	}
	219
de34b538	220	//_____________________________________________________________________________________________
	221	void AliMillePedeRecord::ExpandDtBuffer(Int_t bfsize)
	222	{
	223	// add extra space for derivatives data
	224	bfsize = TMath::Max(bfsize,GetDtBufferSize());
	225	Int_t *tmpI = new Int_t[bfsize];
	226	memcpy(tmpI,fIndex, fSize*sizeof(Int_t));
6fefa5ce	227	delete [] fIndex;
de34b538	228	fIndex = tmpI;
	229	//
	230	Double_t *tmpD = new Double_t[bfsize];
	231	memcpy(tmpD,fValue, fSize*sizeof(Double_t));
6fefa5ce	232	delete [] fValue;
de34b538	233	fValue = tmpD;
	234	//
	235	SetDtBufferSize(bfsize);
	236	}
	237
	238	//_____________________________________________________________________________________________
	239	void AliMillePedeRecord::ExpandGrBuffer(Int_t bfsize)
	240	{
	241	// add extra space for groupID data
	242	bfsize = TMath::Max(bfsize,GetGrBufferSize());
a8c5b94c	243	UShort_t *tmpI = new UShort_t[bfsize];
a8c5b94c	244	memcpy(tmpI,fGroupID, fNGroups*sizeof(UShort_t));
6fefa5ce	245	delete [] fGroupID;
de34b538	246	fGroupID = tmpI;
a8c5b94c	247	for (int i=fNGroups;i<bfsize;i++) fGroupID[i] = 0;
de34b538	248	//
	249	SetGrBufferSize(bfsize);
	250	}
	251
	252	//_____________________________________________________________________________________________
	253	void AliMillePedeRecord::MarkGroup(Int_t id)
	254	{
	255	// mark the presence of the detector group
a8c5b94c	256	id++; // groupID is stored as realID+1
de34b538	257	if (fNGroups>0 && fGroupID[fNGroups-1]==id) return; // already there
	258	if (fNGroups>=GetGrBufferSize()) ExpandGrBuffer(2*(fNGroups+1));
	259	fGroupID[fNGroups++] = id;
	260	}
	261