STEER/STEER/AliMillePedeRecord.cxx

   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() :
  27 fSize(0),fNGroups(0),fRunID(0),fGroupID(0),fIndex(0),fValue(0),fWeight(1) {SetUniqueID(0);}
  28
  29 //_____________________________________________________________________________________________
  30 AliMillePedeRecord::AliMillePedeRecord(const AliMillePedeRecord& src) :
  31   TObject(src),fSize(src.fSize),fNGroups(src.fNGroups),fRunID(src.fRunID),fGroupID(0),fIndex(0),fValue(0),fWeight(src.fWeight)
  32 {
  33   // copy ct-r
  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));
  38   fGroupID = new UShort_t[GetGrBufferSize()];
  39   memcpy(fGroupID,src.fGroupID,GetGrBufferSize()*sizeof(UShort_t));
  40 }
  41
  42 //_____________________________________________________________________________________________
  43 AliMillePedeRecord& AliMillePedeRecord::operator=(const AliMillePedeRecord& rhs)
  44 {
  45   // assignment op-r
  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     }
  54     fWeight = rhs.fWeight;
  55     fRunID = rhs.fRunID;
  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 {
  67   // reset all
  68   fSize = 0;
  69   for (int i=fNGroups;i--;) fGroupID[i] = 0;
  70   fNGroups = 0;
  71   fRunID = 0;
  72   fWeight = 1.;
  73 }
  74
  75 //_____________________________________________________________________________________________
  76 void AliMillePedeRecord::Print(const Option_t *) const
  77 {
  78   // print itself
  79   if (!fSize) {AliInfo("No data"); return;}
  80   int cnt=0,point=0;
  81   //
  82   if (fNGroups) printf("Groups: ");
  83   for (int i=0;i<fNGroups;i++) printf("%4d |",GetGroupID(i));
  84   printf("Run: %9d Weight: %+.2e\n",fRunID,fWeight);
  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
 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 += resid*weight*derGlo[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
 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     //
 209     cnt++;
 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
 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));
 227   delete [] fIndex;
 228   fIndex = tmpI;
 229   //
 230   Double_t *tmpD = new Double_t[bfsize];
 231   memcpy(tmpD,fValue, fSize*sizeof(Double_t));
 232   delete [] fValue;
 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());
 243   UShort_t *tmpI = new UShort_t[bfsize];
 244   memcpy(tmpI,fGroupID, fNGroups*sizeof(UShort_t));
 245   delete [] fGroupID;
 246   fGroupID = tmpI;
 247   for (int i=fNGroups;i<bfsize;i++) fGroupID[i] = 0;
 248   //
 249   SetGrBufferSize(bfsize);
 250 }
 251
 252 //_____________________________________________________________________________________________
 253 void AliMillePedeRecord::MarkGroup(Int_t id)
 254 {
 255   // mark the presence of the detector group
 256   id++; // groupID is stored as realID+1
 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