CTP raw data simulation changed to Run2 payload.

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

diff --git a/STEER/STEER/AliMillePede2.cxx b/STEER/STEER/AliMillePede2.cxx
index 1c9a1fb10693e186b96be3403ca98e7b0653f6d1..c47346710686261ffaa8bc425ab075b36447ae49 100644 (file)
--- a/STEER/STEER/AliMillePede2.cxx
+++ b/STEER/STEER/AliMillePede2.cxx
@@ -31,9 +31,13 @@
 #include <fcntl.h>
 #include <fstream>
 
-using namespace std;
+//#define _DUMP_EQ_BEFORE_
+//#define _DUMP_EQ_AFTER_
 
+//#define _DUMPEQ_BEFORE_
+//#define _DUMPEQ_AFTER_ 
 
+using std::ifstream;
 ClassImp(AliMillePede2)
 
 Bool_t   AliMillePede2::fgInvChol        = kTRUE;     // Invert global matrix with Cholesky solver
@@ -108,6 +112,8 @@ AliMillePede2::AliMillePede2()
   fCurrRecDataID(0),
   fCurrRecConstrID(0),
   fLocFitAdd(kTRUE),
+  fUseRecordWeight(kTRUE),
+  fMinRecordLength(1),
   fSelFirst(1),
   fSelLast(-1),
   fRejRunList(0),
@@ -115,7 +121,9 @@ AliMillePede2::AliMillePede2()
   fAccRunListWgh(0),
   fRunWgh(1),
   fkReGroup(0)
-{}
+{
+  fWghScl[0] = fWghScl[1] = -1;
+}
 
 //_____________________________________________________________________________________________
 AliMillePede2::AliMillePede2(const AliMillePede2& src) : 
@@ -130,6 +138,8 @@ AliMillePede2::AliMillePede2(const AliMillePede2& src) :
   fDataRecFName(0),fRecord(0),fDataRecFile(0),
   fTreeData(0),fRecFileStatus(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
   fCurrRecConstrID(0),fLocFitAdd(kTRUE),
+  fUseRecordWeight(kTRUE),
+  fMinRecordLength(1),
   fSelFirst(1),
   fSelLast(-1),
   fRejRunList(0),
@@ -137,7 +147,11 @@ AliMillePede2::AliMillePede2(const AliMillePede2& src) :
   fAccRunListWgh(0),
   fRunWgh(1),
   fkReGroup(0)
-{printf("Dummy\n");}
+{
+  fWghScl[0] = src.fWghScl[0];
+  fWghScl[1] = src.fWghScl[1];
+  printf("Dummy\n");
+}
 
 //_____________________________________________________________________________________________
 AliMillePede2::~AliMillePede2() 
@@ -230,6 +244,8 @@ Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut,
     fParamGrID[i] = -1;
   }
   //
+  fWghScl[0] = -1; 
+  fWghScl[1] = -1;
   return 1;
 }
 
@@ -543,14 +559,19 @@ Int_t AliMillePede2::LocalFit(double *localParams)
     //
     nPoints++;
   }
+  if (fMinRecordLength>0 && nPoints < fMinRecordLength) return 0; // ignore
+  //
   double vl;
   //
-  double gloWgh = fRecord->GetWeight()*fRunWgh; // global weight for this set
+  double gloWgh = fRunWgh;
+  if (fUseRecordWeight) 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 )*gloWgh;    
+    double  weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+    int odd = (ip&0x1);
+    if (fWghScl[odd]>0) weight *= fWghScl[odd];
     double *derLoc = fRecord->GetValue()+refLoc[ip];
     double *derGlo = fRecord->GetValue()+refGlo[ip];
     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
@@ -610,7 +631,9 @@ 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 )*gloWgh;    
+    double  weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+    int odd = (ip&0x1);
+    if (fWghScl[odd]>0) weight *= fWghScl[odd];
     double *derLoc = fRecord->GetValue()+refLoc[ip];
     double *derGlo = fRecord->GetValue()+refGlo[ip];
     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
@@ -665,7 +688,9 @@ 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 )*gloWgh;    
+    double  weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+    int odd = (ip&0x1);
+    if (fWghScl[odd]>0) weight *= fWghScl[odd];
     double *derLoc = fRecord->GetValue()+refLoc[ip];
     double *derGlo = fRecord->GetValue()+refGlo[ip];
     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
@@ -1049,12 +1074,60 @@ Int_t AliMillePede2::GlobalFitIteration()
   //
   sws.Start();
 
+#ifdef _DUMP_EQ_BEFORE_
+  const char* faildumpB = Form("mp2eq_before%d.dat",fIter);
+  int defoutB = dup(1);
+  if (defoutB<0) {AliFatal("Failed on dup"); exit(1);}
+  int slvDumpB = open(faildumpB, O_RDWR|O_CREAT, 0666);
+  if (slvDumpB>=0) {
+    dup2(slvDumpB,1);
+    printf("Solving%d for %d params\n",fIter,fNGloSize);
+    matCGlo.Print("10");
+    for (int i=0;i<fNGloSize;i++) printf("b%2d : %+.10f\n",i,fVecBGlo[i]);
+  }
+  dup2(defoutB,1);
+  close(slvDumpB);
+  close(defoutB);
+
+#endif
   /*
   printf("Solving:\n");
   matCGlo.Print("l");
   for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e\n",i,fVecBGlo[i]);
   */
+#ifdef _DUMPEQ_BEFORE_  
+  const char* faildumpB = Form("mp2eq_before%d.dat",fIter);
+  int defoutB = dup(1);
+  int slvDumpB = open(faildumpB, O_RDWR|O_CREAT, 0666);
+  dup2(slvDumpB,1);
+  //
+  printf("#Equation before step %d\n",fIter);
+  fMatCGlo->Print("10");
+  printf("#RHS/STAT : NGlo:%d NGloSize:%d\n",fNGloPar,fNGloSize);
+  for (int i=0;i<fNGloSize;i++) printf("%d %+.10f %d\n",i,fVecBGlo[i],fProcPnt[i]);
+  //
+  dup2(defoutB,1);
+  close(slvDumpB);
+  close(defoutB);
+#endif
+  //
   fGloSolveStatus = SolveGlobalMatEq();                     // obtain solution for this step
+#ifdef _DUMPEQ_AFTER_  
+  const char* faildumpA = Form("mp2eq_after%d.dat",fIter);
+  int defoutA = dup(1);
+  int slvDumpA = open(faildumpA, O_RDWR|O_CREAT, 0666);
+  dup2(slvDumpA,1);
+  //
+  printf("#Matrix after step %d\n",fIter);
+  fMatCGlo->Print("10");
+  printf("#RHS/STAT : NGlo:%d NGloSize:%d\n",fNGloPar,fNGloSize);
+  for (int i=0;i<fNGloSize;i++) printf("%d %+.10f %d\n",i,fVecBGlo[i],fProcPnt[i]);
+  //
+  dup2(defoutA,1);
+  close(slvDumpA);
+  close(defoutA);
+#endif
+  //
   sws.Stop();
   printf("Solve %d |",fIter); sws.Print();
   //
@@ -1063,6 +1136,22 @@ Int_t AliMillePede2::GlobalFitIteration()
   if (fGloSolveStatus==kFailed) return 0;
   //
   for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
+
+#ifdef _DUMP_EQ_AFTER_
+  const char* faildumpA = Form("mp2eq_after%d.dat",fIter);
+  int defoutA = dup(1);
+  if (defoutA<0) {AliFatal("Failed on dup"); exit(1);}
+  int slvDumpA = open(faildumpA, O_RDWR|O_CREAT, 0666);
+  if (slvDumpA>=0) {
+    dup2(slvDumpA,1);
+    printf("Solving%d for %d params\n",fIter,fNGloSize);
+    matCGlo.Print("10");
+    for (int i=0;i<fNGloSize;i++) printf("b%2d : %+.10f\n",i,fVecBGlo[i]);
+  }
+  dup2(defoutA,1);
+  close(slvDumpA);
+  close(defoutA);
+#endif
   //
   /*
   printf("Solved:\n");
@@ -1070,7 +1159,7 @@ Int_t AliMillePede2::GlobalFitIteration()
   for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e (->%+e)\n",i,fVecBGlo[i], fDeltaPar[i]);
   */
 
-  //  PrintGlobalParameters();
+  PrintGlobalParameters();
   return 1;
 }
 
@@ -1131,6 +1220,7 @@ Int_t AliMillePede2::SolveGlobalMatEq()
       //
       dup2(defout,1);
       close(slvDump);
+      close(defout);
       printf("#Dumped failed matrix and RHS to %s\n",faildump);
     }
     else AliInfo("Failed on file open for matrix dumping");
@@ -1138,6 +1228,7 @@ Int_t AliMillePede2::SolveGlobalMatEq()
     return kFailed;
   }
   for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i];
+  //
   return kNoInversion;
   //
 }
@@ -1274,13 +1365,19 @@ Bool_t AliMillePede2::IsRecordAcceptable()
       prevAns = kTRUE;
       for (int i=n;i--;) if (runID == (*fRejRunList)[i]) {
          prevAns = kFALSE; 
-         printf("New Run to reject: %ld -> %d\n",runID,prevAns);
+         AliInfo(Form("New Run to reject: %ld",runID));
          break;
        }
     }
     else if (fAccRunList && (n=fAccRunList->GetSize())) {     // is run specifically selected
       prevAns = kFALSE;
-      for (int i=n;i--;) if (runID == (*fAccRunList)[i]) {prevAns = kTRUE; fRunWgh = (*fAccRunListWgh)[i]; break;}
+      for (int i=n;i--;) if (runID == (*fAccRunList)[i]) {
+        prevAns = kTRUE; 
+       if (fAccRunListWgh) fRunWgh = (*fAccRunListWgh)[i]; 
+        AliInfo(Form("New Run to accept explicitly: %ld, weight=%f",runID,fRunWgh));
+       break;
+      }
+      if (!prevAns) AliInfo(Form("New Run is not in the list to accept: %ld",runID)); 
     }
   }
   //