doubles in the AliMillePedeRecord changed to double32, many changes in

MP2 related to treatment of constraints and low-statistics parameters

diff --git a/STEER/AliMillePede2.cxx b/STEER/AliMillePede2.cxx
index a18dfca5afa9e833428cf15703e714b9f30b857b..4963e0d17023c26a3c0c170899cbd6568e205888 100644 (file)
--- a/STEER/AliMillePede2.cxx
+++ b/STEER/AliMillePede2.cxx
@@ -82,6 +82,7 @@ AliMillePede2::AliMillePede2()
   fRecord(0),
   fDataRecFile(0),  
   fTreeData(0),
+  fRecFileStatus(0),
   //
   fConstrRecFName("/tmp/mp2_constraints_records.root"),
   fTreeConstr(0),
@@ -101,7 +102,7 @@ AliMillePede2::AliMillePede2(const AliMillePede2& src) :
   fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0),
   fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0),
   fDataRecFName(0),fRecord(0),fDataRecFile(0),
-  fTreeData(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
+  fTreeData(0),fRecFileStatus(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
   fCurrRecConstrID(0),fLocFitAdd(kTRUE)
 {printf("Dummy\n");}
 
@@ -215,12 +216,12 @@ Bool_t AliMillePede2::InitDataRecStorage(Bool_t read)
   if (!fRecord) fRecord = new AliMillePedeRecord();
   //
   fDataRecFile = TFile::Open(GetDataRecFName(),read ? "":"recreate");
-  if (!fDataRecFile) {AliInfo(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;}
+  if (!fDataRecFile) {AliFatal(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;}
   //
   AliInfo(Form("File %s used for derivatives records",GetDataRecFName()));
   if (read) {
     fTreeData = (TTree*)fDataRecFile->Get("AliMillePedeRecords_Data");
-    if (!fTreeData) {AliInfo(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;}
+    if (!fTreeData) {AliFatal(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;}
     fTreeData->SetBranchAddress("Record_Data",&fRecord);
     AliInfo(Form("Found %d derivatives records",fTreeData->GetEntries()));
   }
@@ -229,6 +230,7 @@ Bool_t AliMillePede2::InitDataRecStorage(Bool_t read)
     fTreeData->Branch("Record_Data","AliMillePedeRecord",&fRecord,32000,99);
   }
   fCurrRecDataID = -1;
+  fRecFileStatus = read ? 1:2;
   //
   return kTRUE;
 }
@@ -277,6 +279,7 @@ void AliMillePede2::CloseDataRecStorage()
     delete fDataRecFile;
     fDataRecFile = 0;
   }
+  fRecFileStatus = 0;
   //
 }
 
@@ -315,10 +318,17 @@ Bool_t AliMillePede2::ReadNextRecordConstraint()
   return kTRUE;
 }
 
+//_____________________________________________________________________________________________
+void AliMillePede2::SetRecordWeight(double wgh)
+{
+  if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
+  fRecord->SetWeight(wgh);
+}
+
 //_____________________________________________________________________________________________
 void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma)
 {
-  if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
+  if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
   //
   // write data of single measurement
   if (lSigma<=0.0) { // If parameter is fixed, then no equation
@@ -353,7 +363,7 @@ void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *in
     return;
   }
   //
-  if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
+  if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
   //
   fRecord->AddResidual(lMeas);
   //
@@ -450,17 +460,18 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   }
   double vl;
   //
+  double gloWgh = fRecord->GetWeight(); // global weight for this set
   Int_t maxLocUsed = 0;
   //
   for (int ip=nPoints;ip--;) {  // Transfer the measurement records to matrices
     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
-    double  weight = fRecord->GetValue( refGlo[ip]-1 );    
+    double  weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;    
     double *derLoc = fRecord->GetValue()+refLoc[ip];
     double *derGlo = fRecord->GetValue()+refGlo[ip];
     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
     int    *indGlo = fRecord->GetIndex()+refGlo[ip];
     //
-    for (int i=nrefGlo[ip];i--;) {       // suppress the global part (only relevant with iterations)
+    for (int i=nrefGlo[ip];i--;) {      // suppress the global part (only relevant with iterations)
       int iID = indGlo[i];              // Global param indice
       if (fSigmaPar[iID]<=0.) continue;                                    // fixed parameter RRRCheck
       if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);  // linear parameter
@@ -480,7 +491,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   //
   // first try to solve by faster Cholesky decomposition, then by Gaussian elimination
   if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) {
-    AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
+    AliInfo("Failed to solve locals by Cholesky, trying Gaussian Elimination");
     if (!matCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { 
       AliInfo("Failed to solve locals by Gaussian Elimination, skip...");
       matCLoc.Print("d");
@@ -500,7 +511,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   //
   for (int ip=nPoints;ip--;) {  // Calculate residuals
     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
-    double  weight = fRecord->GetValue( refGlo[ip]-1 );    
+    double  weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;    
     double *derLoc = fRecord->GetValue()+refLoc[ip];
     double *derGlo = fRecord->GetValue()+refGlo[ip];
     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
@@ -529,6 +540,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
     nEq++;                        // number of equations                       
   } // end of Calculate residuals
   //
+  lChi2 /= gloWgh;  
   int nDoF = nEq-maxLocUsed;
   lChi2 = (nDoF>0) ? lChi2/nDoF : 0;  // Chi^2/dof  
   //
@@ -554,7 +566,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   //
   for (int ip=nPoints;ip--;) {  // Update matrices
     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
-    double  weight = fRecord->GetValue( refGlo[ip]-1 );    
+    double  weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;    
     double *derLoc = fRecord->GetValue()+refLoc[ip];
     double *derGlo = fRecord->GetValue()+refGlo[ip];
     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
@@ -832,7 +844,7 @@ Int_t AliMillePede2::GlobalFitIteration()
     if (fProcPnt[i]<1) {
       fNGloFix++; 
       fVecBGlo[i] = 0.;
-      matCGlo.DiagElem(i) = float(fNLocEquations*fNLocEquations);
+      matCGlo.DiagElem(i) = 1.;//float(fNLocEquations*fNLocEquations);
     }
     else matCGlo.DiagElem(i) += (fgWeightSigma ? fProcPnt[i] : 1.)/(fSigmaPar[i]*fSigmaPar[i]);
   }

diff --git a/STEER/AliMillePede2.h b/STEER/AliMillePede2.h
index e949e24d21c1e988754ef0b08e1efe6dd064b36f..606e0d9fecc494d15734baeeb75cbb6e536ccd05 100644 (file)
--- a/STEER/AliMillePede2.h
+++ b/STEER/AliMillePede2.h
@@ -123,6 +123,7 @@ class AliMillePede2: public TObject
   void                 SetGlobalConstraint(const int *indgb,double *dergb, int ngb, double val, double sigma=0);
   //
   // processing of the local measurement
+  void                 SetRecordWeight(double wgh);
   void                 SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma);
   void                 SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc, 
                                        double *derlc,int nlc,double lMeas,double lSigma);
@@ -210,6 +211,7 @@ class AliMillePede2: public TObject
   AliMillePedeRecord   *fRecord;                         // Buffer of measurements records
   TFile                *fDataRecFile;                    // File of processed measurements records
   TTree                *fTreeData;                       // Tree of processed measurements records
+  Int_t                 fRecFileStatus;                  // state of the record file (0-no, 1-read, 2-rw)
   //
   TString               fConstrRecFName;                 // Name of File for constraints records               
   TTree                *fTreeConstr;                     // Tree of constraint records

diff --git a/STEER/AliMillePedeRecord.cxx b/STEER/AliMillePedeRecord.cxx
index 4cba214c5ff3e1c45da0b61ebdb7f2051220b77a..c6ff01ff4f17cb190b2989cca550a0412c0a9886 100644 (file)
--- a/STEER/AliMillePedeRecord.cxx
+++ b/STEER/AliMillePedeRecord.cxx
@@ -28,14 +28,14 @@ fSize(0),fNGroups(0),fGroupID(0),fIndex(0),fValue(0) {SetUniqueID(0);}
 
 //_____________________________________________________________________________________________
 AliMillePedeRecord::AliMillePedeRecord(const AliMillePedeRecord& src) : 
-  TObject(src),fSize(src.fSize),fNGroups(src.fNGroups),fGroupID(0),fIndex(0),fValue(0)
+  TObject(src),fSize(src.fSize),fNGroups(src.fNGroups),fGroupID(0),fIndex(0),fValue(0),fWeight(src.fWeight)
 {
   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 Int_t[GetGrBufferSize()];
-  memcpy(fGroupID,src.fGroupID,GetGrBufferSize()*sizeof(Int_t));
+  fGroupID = new UShort_t[GetGrBufferSize()];
+  memcpy(fGroupID,src.fGroupID,GetGrBufferSize()*sizeof(UShort_t));
 }
 
 //_____________________________________________________________________________________________
@@ -49,6 +49,7 @@ AliMillePedeRecord& AliMillePedeRecord::operator=(const AliMillePedeRecord& rhs)
       rhs.GetIndexValue(i,ind,val);
       AddIndexValue(ind,val);
     }
+    fWeight = rhs.fWeight;
     for (int i=0;i<rhs.GetNGroups();i++) MarkGroup(rhs.GetGroupID(i));
   }
   return *this;
@@ -61,8 +62,9 @@ AliMillePedeRecord::~AliMillePedeRecord() {delete[] fIndex; delete[] fValue; del
 void AliMillePedeRecord::Reset()
 {
   fSize = 0;
-  for (int i=fNGroups;i--;) fGroupID[i] = -1;
+  for (int i=fNGroups;i--;) fGroupID[i] = 0;
   fNGroups = 0;
+  fWeight = 1.;
 }
 
 //_____________________________________________________________________________________________
@@ -72,7 +74,8 @@ void AliMillePedeRecord::Print(const Option_t *) const
   int cnt=0,point=0;
   //  
   if (fNGroups) printf("Groups: ");
-  for (int i=0;i<fNGroups;i++) printf("%4d |",GetGroupID(i)); printf("\n");
+  for (int i=0;i<fNGroups;i++) printf("%4d |",GetGroupID(i)); 
+  printf(" Weight: %+.2e\n",fWeight);
   while(cnt<fSize) {
     //
     Double_t resid = fValue[cnt++];
@@ -119,11 +122,11 @@ void AliMillePedeRecord::ExpandGrBuffer(Int_t bfsize)
 {
   // add extra space for groupID data 
   bfsize = TMath::Max(bfsize,GetGrBufferSize());
-  Int_t *tmpI = new Int_t[bfsize];
-  memcpy(tmpI,fGroupID, fNGroups*sizeof(Int_t));
-  delete fGroupID;
+  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] = -1;
+  for (int i=fNGroups;i<bfsize;i++) fGroupID[i] = 0;
   //
   SetGrBufferSize(bfsize);
 }
@@ -132,14 +135,9 @@ void AliMillePedeRecord::ExpandGrBuffer(Int_t 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;  
 }
 
-//_____________________________________________________________________________________________
-Bool_t AliMillePedeRecord::IsGroupPresent(Int_t id) const
-{
-  for (int i=fNGroups;i--;) if (GetGroupID(i)==id) return kTRUE;
-  return kFALSE;
-}

diff --git a/STEER/AliMillePedeRecord.h b/STEER/AliMillePedeRecord.h
index 2dee67f516d2c1d37c03c732cc92ec07e4d44c25..94e8e073cc01f155070bf41af50b12e8ce841a96 100644 (file)
--- a/STEER/AliMillePedeRecord.h
+++ b/STEER/AliMillePedeRecord.h
@@ -39,15 +39,17 @@ class AliMillePedeRecord : public TObject
   void       AddIndexValue(Int_t ind, Double_t val);
   void       AddResidual(Double_t val)                             {AddIndexValue(-1,val);}
   void       AddWeight(Double_t val)                               {AddIndexValue(-2,val);}
+  void       SetWeight(Double_t w=1)                               {fWeight = w;}
   Bool_t     IsResidual(Int_t i)                             const {return fIndex[i]==-1;}
   Bool_t     IsWeight(Int_t i)                               const {return fIndex[i]==-2;}
   //
   Double_t  *GetValue()                                      const {return fValue;}
   Double_t   GetValue(Int_t i)                               const {return fValue[i];}
+  Double_t   GetWeight()                                     const {return fWeight;}
   //
   void       MarkGroup(Int_t id);
   Int_t      GetNGroups()                                    const {return fNGroups;}
-  Int_t      GetGroupID(Int_t i)                             const {return fGroupID[i];}
+  Int_t      GetGroupID(Int_t i)                             const {return fGroupID[i]-1;}
   Bool_t     IsGroupPresent(Int_t id)                        const;
   //
  protected:
@@ -61,18 +63,30 @@ class AliMillePedeRecord : public TObject
  protected:
   Int_t      fSize;                             // size of the record
   Int_t      fNGroups;                          // number of groups (e.g. detectors) contributing
-  Int_t   *  fGroupID;                          //[fNGroups] groups id's (in increasing order)
+  UShort_t*  fGroupID;                          //[fNGroups] groups id's+1 (in increasing order)
   Int_t   *  fIndex;                            //[fSize] index of variables
-  Double_t*  fValue;                            //[fSize] array of values: derivs,residuals
+  Double32_t* fValue;                           //[fSize] array of values: derivs,residuals
+  Double32_t  fWeight;                          //global weight for the record
   //
-  ClassDef(AliMillePedeRecord,1)                // Record of track residuals and local/global deriavtives
+  ClassDef(AliMillePedeRecord,2)                // Record of track residuals and local/global deriavtives
 };
 
+//_____________________________________________________________________________________________
 inline void  AliMillePedeRecord::AddIndexValue(Int_t ind, Double_t val) 
 {
+  // add new pair of index/value
   if (fSize>=GetDtBufferSize()) ExpandDtBuffer(2*(fSize+1));
   fIndex[fSize]=ind; 
   fValue[fSize++]=val;
 }
 
+//_____________________________________________________________________________________________
+inline Bool_t AliMillePedeRecord::IsGroupPresent(Int_t id) const
+{
+  // check if group is defined
+  id++;
+  for (int i=fNGroups;i--;) if (fGroupID[i]==id) return kTRUE;
+  return kFALSE;
+}
+
 #endif