o add SC scaling factor calculation and usage in
authorwiechula <wiechula@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 Oct 2013 08:26:00 +0000 (08:26 +0000)
committerwiechula <wiechula@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 Oct 2013 08:26:00 +0000 (08:26 +0000)
  - the correction class
  - the generator
  - the reconstructor

TPC/Base/AliTPCCorrection.h
TPC/Base/AliTPCCorrectionLookupTable.cxx
TPC/Base/AliTPCCorrectionLookupTable.h
TPC/Upgrade/AliToyMCEvent.cxx
TPC/Upgrade/AliToyMCEvent.h
TPC/Upgrade/AliToyMCEventGenerator.cxx
TPC/Upgrade/AliToyMCEventGenerator.h
TPC/Upgrade/AliToyMCReconstruction.cxx

index 930643b04d0f2b5d09c31b7eebb65c3bac43b2d6..348a2d8a951cadb3a7745080ae79dc1741c72c13 100644 (file)
@@ -54,6 +54,10 @@ public:
   virtual void Init();
   virtual void Update(const TTimeStamp &timeStamp);
 
+  // map scaling
+  virtual void    SetCorrScaleFactor(Float_t /*val*/) { ; }
+  virtual Float_t GetCorrScaleFactor() const { return 1.; }
+  
   // convenience functions
   virtual void Print(Option_t* option="") const;
  
index ab568650b654c657c678e8ffcca71c57ed766fca..31e394cab009f8511f97a6905b713b957b5e8951 100644 (file)
@@ -39,6 +39,7 @@ AliTPCCorrectionLookupTable::AliTPCCorrectionLookupTable()
 , fNR(0)
 , fNPhi(0)
 , fNZ(0)
+, fCorrScaleFactor(-1)
 , fFillCorrection(kTRUE)
 , fLimitsR()
 , fLimitsPhi()
@@ -79,6 +80,12 @@ void AliTPCCorrectionLookupTable::GetCorrection(const Float_t x[],const Short_t
   // Get interplolated correction
   //
   GetInterpolation(x,roc,dx,fLookUpDxCorr,fLookUpDyCorr,fLookUpDzCorr);
+
+  if (fCorrScaleFactor>0) {
+    dx[0]*=fCorrScaleFactor;
+    dx[1]*=fCorrScaleFactor;
+    dx[2]*=fCorrScaleFactor;
+  }
 }
 
 //_________________________________________________________________________________________
index 4a39a488e6be1edbd0d78f40b189a5f730914dbe..af2d651465e4402e3dbe2a653a5d8c43426ce488 100644 (file)
@@ -37,6 +37,9 @@ public:
   const TVectorD& GetLimitsPhi() const { return fLimitsPhi; }
   const TVectorD& GetLimitsZ()   const { return fLimitsZ; }
   
+  void SetCorrScaleFactor(Float_t    val) { fCorrScaleFactor = val; }
+  Float_t    GetCorrScaleFactor() const { return fCorrScaleFactor; }
+  
 private:
 
   // sizes of lookup tables
@@ -44,6 +47,9 @@ private:
   Int_t     fNR;                   // number of rows (r) used for lookup table
   Int_t     fNPhi;                 // number of phi slices used for lookup table
   Int_t     fNZ;                   // number of columns (z) used for lookup table
+
+  Float_t   fCorrScaleFactor;      // overall scaling factor for the correction
+  
   Bool_t    fFillCorrection;       // whether to also fill the correction tables
   //
   TVectorD  fLimitsR;              // bin limits in row direction
index 3ec158d10d7f33264a25c9a85112cfdc86dd5fe6..6b3705b34f0758b90973e4d596e93b89e7c1350a 100644 (file)
@@ -11,6 +11,8 @@ AliToyMCEvent::AliToyMCEvent()
   ,fX(-1000.)
   ,fY(-1000.)
   ,fZ(-1000.)
+  ,fSCscale(-1.)
+  ,fSCscaleChi2(0)
   ,fTracks("AliToyMCTrack")
 {
   fEventNumber = fgEvCounter;
@@ -26,6 +28,8 @@ AliToyMCEvent::AliToyMCEvent(const AliToyMCEvent &event)
   ,fX(event.fX)
   ,fY(event.fY)
   ,fZ(event.fZ)
+  ,fSCscale(event.fSCscale)
+  ,fSCscaleChi2(event.fSCscaleChi2)
   ,fTracks(event.fTracks)
 {
   //
index dc70ed5a8aed85e146a63ec6e43f21c07f406c4d..c08d5a602574d640f2c4d44887ef4ab4dbdeb548 100644 (file)
@@ -39,22 +39,31 @@ class AliToyMCEvent : public TObject {
 
   void SetEventType(EEventType type) { fEventType=type; }
   EEventType GetEventType() const    { return fEventType; }
+
+  void     SetSCscale(Float_t  val)      { fSCscale = val;      }
+  Float_t  GetSCscale() const            { return fSCscale;     }
+
+  void     SetSCscaleChi2(Float_t  val)  { fSCscaleChi2 = val;  }
+  Float_t  GetSCscaleChi2() const        { return fSCscaleChi2; }
   
  private:
-  static Int_t fgEvCounter;
+  static Int_t fgEvCounter;      // static counter
     
-  UInt_t fEventNumber;
+  UInt_t fEventNumber;           // event number
 
-  EEventType fEventType;
+  EEventType fEventType;         // type of the event
     
-  Float_t fT0;
-  Float_t fX;
-  Float_t fY;
-  Float_t fZ;
+  Float_t fT0;                   // Interaction time of the event
+  Float_t fX;                    // x-vertex position
+  Float_t fY;                    // y-vertex position
+  Float_t fZ;                    // z-vertex position
 
-  TClonesArray fTracks;
+  Float_t fSCscale;              // scaling parameter for space charge correction
+  Float_t fSCscaleChi2;          // chi2 of scaling parameter for space charge correction
+  
+  TClonesArray fTracks;          // array of tracks
 
-  ClassDef(AliToyMCEvent, 1);
+  ClassDef(AliToyMCEvent, 2);    // MC event class
 
 };
 
index e39311b691fb353d6c34269e21750dc57db87e04..c2acbcd0f268c0c9acf346c6bbfb371c78295a3a 100644 (file)
@@ -10,6 +10,7 @@
 #include <TROOT.h>
 #include <TSystem.h>
 #include <TObjArray.h>
+#include <TLinearFitter.h>
 
 #include <AliLog.h>
 #include <AliTPCROC.h>
@@ -38,6 +39,7 @@ AliToyMCEventGenerator::AliToyMCEventGenerator()
   ,fEvent(0x0)
   ,fCurrentTrack(0)
   ,fTPCCorrection(0x0)
+  ,fTPCCorrectionAv(0x0)
   ,fSCList(0x0)
   ,fSCListFile()
   ,fCorrectionFile("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.lookup.root")
@@ -60,6 +62,7 @@ AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen
   ,fEvent(0x0)
   ,fCurrentTrack(0)
   ,fTPCCorrection(gen.fTPCCorrection)
+  ,fTPCCorrectionAv(0x0)
   ,fSCList(gen.fSCList)
   ,fSCListFile(gen.fSCListFile)
   ,fCorrectionFile(gen.fCorrectionFile)
@@ -621,20 +624,7 @@ void AliToyMCEventGenerator::InitSpaceCharge()
   
   AliInfo(Form("Using space charge map file: '%s'",fCorrectionFile.Data()));
 
-  TString corrName("map");
-
-  // allow for specifying an object name for the AliTPCCorrection in the file name
-  // separated by a ':'
-  TObjArray *arr=fCorrectionFile.Tokenize(":");
-  if (arr->GetEntriesFast()>1) {
-    fCorrectionFile=arr->At(0)->GetName();
-    corrName=arr->At(1)->GetName();
-  }
-  delete arr;
-  
-  
-  TFile f(fCorrectionFile.Data());
-  fTPCCorrection=(AliTPCSpaceCharge3D*)f.Get("map");
+  SetCorrectionFromFile(fCorrectionFile, fTPCCorrection);
 
   if (fOutTree){
     AliInfo("Attaching space charge map file name to the tree");
@@ -673,6 +663,10 @@ void AliToyMCEventGenerator::InitSpaceChargeList()
     fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
   }
   
+  // check for an average map
+  if (!fCorrectionFile.IsNull()) {
+    SetCorrectionFromFile(fCorrectionFile, fTPCCorrectionAv);
+  }
   
   // case of non preread
   // store the names of the files
@@ -689,11 +683,10 @@ void AliToyMCEventGenerator::InitSpaceChargeList()
   fSCList->Delete();
   
   for (Int_t ifile=0; ifile<arr->GetEntriesFast(); ++ifile) {
-    TFile f(arr->At(ifile)->GetName());
-    if (!f.IsOpen()||f.IsZombie()) continue;
-    AliTPCSpaceCharge3D *sc=(AliTPCSpaceCharge3D*)f.Get("map");
+    TString scname=arr->At(ifile)->GetName();
+    AliTPCCorrection *sc=0x0;
+    SetCorrectionFromFile(scname, sc);
     if (!sc) continue;
-    sc->SetName(f.GetName());
     fSCList->Add(sc);
   }
 }
@@ -725,16 +718,100 @@ void AliToyMCEventGenerator::IterateSC(Int_t ipos)
   delete fTPCCorrection;
   fTPCCorrection=0x0;
 
-  TFile f(file);
-  if (!f.IsOpen()||f.IsZombie()) {
-    AliFatal(Form("Could not open SC file '%s'",file.Data()));
+  SetCorrectionFromFile(file, fTPCCorrection);
+  
+  if (!fTPCCorrection) {
+    AliFatal(Form("Could not read SC map from SC file '%s'",file.Data()));
     return;
   }
 
-  fTPCCorrection=(AliTPCSpaceCharge3D*)f.Get("map");
-  if (!fTPCCorrection) {
-    AliFatal(Form("Could not read SC map from SC file '%s'",file.Data()));
+  SetSCScalingFactor();
+}
+
+//________________________________________________________________
+Float_t AliToyMCEventGenerator::GetSCScalingFactor(AliTPCCorrection *corr, AliTPCCorrection *averageCorr, Float_t &chi2)
+{
+  //
+  // calculate the scaling factor
+  // between the fluctuation map and the average map
+  //
+
+  TLinearFitter fitter(2,"pol1");
+  
+  for (Float_t iz=-245; iz<=245; iz+=10) {
+    Short_t roc=(iz>=0)?0:18;
+    for (Float_t ir=86; ir<250; ir+=10) {
+      for (Float_t iphi=0; iphi<TMath::TwoPi(); iphi+=10*TMath::DegToRad()){
+        Float_t x=ir*(Float_t)TMath::Cos(iphi);
+        Float_t y=ir*(Float_t)TMath::Sin(iphi);
+        Float_t x3[3]    = {x,y,iz};
+        Float_t dx3[3]   = {0.,0.,0.};
+        Float_t dx3av[3] = {0.,0.,0.};
+        
+        corr->GetDistortion(x3,roc,dx3);
+        averageCorr->GetDistortion(x3,roc,dx3av);
+
+        Double_t dr   = TMath::Sqrt(dx3[0]*dx3[0]     +  dx3[1]*dx3[1]);
+        Double_t drav = TMath::Sqrt(dx3av[0]*dx3av[0] +  dx3av[1]*dx3av[1]);
+
+        fitter.AddPoint(&drav,dr);
+      }
+    }
+  }
+  
+  fitter.Eval();
+  chi2 = fitter.GetChisquare()/fitter.GetNpoints();
+  return fitter.GetParameter(1);
+}
+
+//________________________________________________________________
+void AliToyMCEventGenerator::SetSCScalingFactor()
+{
+  //
+  // if running with a SC list calculate the scaling factor
+  // between the fluctuation map and the average map
+  //
+
+  if ( !(fTPCCorrection && fTPCCorrectionAv && fEvent) ) return;
+
+  // loop over several z, r and phi bins and find a factor that minimises
+  // the distortions between the current map and the average map
+
+  Float_t chi2   = 0.;
+  Float_t factor = GetSCScalingFactor(fTPCCorrection, fTPCCorrectionAv, chi2);
+
+  fEvent->SetSCscale(factor);
+  fEvent->SetSCscaleChi2(chi2);
+}
+
+//________________________________________________________________
+void AliToyMCEventGenerator::SetCorrectionFromFile(const TString& file, AliTPCCorrection* &corr)
+{
+  //
+  // read the correction from file and set it to corr
+  //
+
+  corr=0x0;
+  TString corrName("map");
+  
+  // allow for specifying an object name for the AliTPCCorrection in the file name
+  // separated by a ':'
+  TObjArray *arr=fCorrectionFile.Tokenize(":");
+  if (arr->GetEntriesFast()>1) {
+    fCorrectionFile=arr->At(0)->GetName();
+    corrName=arr->At(1)->GetName();
+  }
+  delete arr;
+  
+  
+  TFile f(file.Data());
+  if (!f.IsOpen()||f.IsZombie()) {
+    AliError(Form("Could not open SC file '%s'",file.Data()));
     return;
   }
   
+  corr=(AliTPCSpaceCharge3D*)f.Get(corrName.Data());
+  if (corr) corr->SetName(f.GetName());
 }
+
+
index 043d24e888b612a0df8a8b39aa5ea71eb2a3a98c..67d2288ea35f170ace7c8fe6dab16d482c3f357b 100644 (file)
@@ -81,6 +81,8 @@ class AliToyMCEventGenerator : public TObject {
   void   SetPrereadSCList(Bool_t b)      { fPrereadSCList=b;              }
   Bool_t GetPrereadSCList() const        { return fPrereadSCList;         }
   Bool_t HasSCList() const               { return  !fSCListFile.IsNull(); }
+
+  static Float_t GetSCScalingFactor(AliTPCCorrection *corr, AliTPCCorrection *averageCorr, Float_t &chi2);
   
  protected:
   AliTPCParam *fTPCParam;                //! TPC params
@@ -98,7 +100,8 @@ class AliToyMCEventGenerator : public TObject {
   AliToyMCEventGenerator& operator= (const AliToyMCEventGenerator& );
    
   AliTPCCorrection *fTPCCorrection;      //! distortion correction
-
+  AliTPCCorrection *fTPCCorrectionAv;    //! average distortion correction
+  
   TObjArray   *fSCList;                  //! list with
   TString fSCListFile;                   // file with a list of space charge files
   TString fCorrectionFile;               // name of a sinfle SC file
@@ -111,7 +114,9 @@ class AliToyMCEventGenerator : public TObject {
   Bool_t fIsLaser;                       // is a laser event?
   Bool_t fPrereadSCList;                 // preread all SC files from the SC list
 
+  void SetCorrectionFromFile(const TString& file, AliTPCCorrection* &corr);
   void InitSpaceChargeList();
+  void SetSCScalingFactor();
   
   ClassDef(AliToyMCEventGenerator, 1)
   
index 34abc2b67d3d3092ce0b99f1e34a942a7d335c1a..11527834a641f90be91c387d1b1afe03c3f2e439 100644 (file)
@@ -159,6 +159,9 @@ void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
     
     Float_t z0=fEvent->GetZ();
     Float_t t0=fEvent->GetT0();
+
+    // set SC scaling factor
+    fTPCCorrection->SetCorrScaleFactor(fEvent->GetSCscale());
     
     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 //       printf(" > ======  Processing Track %6d ========  \n",itr);