Truncated mean PID method by Xianguo
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 Mar 2012 11:21:32 +0000 (11:21 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 Mar 2012 11:21:32 +0000 (11:21 +0000)
TRD/AliTRDCalibTask.cxx
TRD/AliTRDPreprocessorOffline.cxx
TRD/AliTRDPreprocessorOffline.h
TRD/AliTRDcalibDB.cxx
TRD/AliTRDcalibDB.h
TRD/AliTRDtrackV1.cxx
TRD/AliTRDtrackV1.h
TRD/CMakelibTRDrec.pkg
TRD/TRDrecLinkDef.h

index 52aa79f..6df00b0 100644 (file)
@@ -78,6 +78,7 @@ using namespace std;
 #include "AliLog.h"
 
 #include "AliTRDCalibChamberStatus.h"
+#include "AliTRDdEdxUtils.h"
 
 #include "AliTRDCalibTask.h"
 
@@ -220,6 +221,7 @@ AliTRDCalibTask::~AliTRDCalibTask()
   if(fCH2dTest) delete fCH2dTest;
   if(fPH2dTest) delete fPH2dTest;
   if(fLinearVdriftTest) delete fLinearVdriftTest;
+  AliTRDdEdxUtils::DeleteCalibHist();
   if(fCalDetGain) delete fCalDetGain;
   
   if(fSelectedTrigger) {
@@ -340,6 +342,8 @@ void AliTRDCalibTask::UserCreateOutputObjects()
   fAbsoluteGain->Sumw2();
   fListHist->Add(fAbsoluteGain);
   
+  AliTRDdEdxUtils::PrintControl();
+  AliTRDdEdxUtils::IniCalibHist(fListHist, kTRUE);
   /////////////////////////////////////////
   // First debug level
   ///////////////////////////////////////
@@ -561,7 +565,7 @@ void AliTRDCalibTask::UserExec(Option_t *)
   //cout << "event = " << fCounter << endl;
   
   //printf("Counter %d\n",fCounter);
-  
+
   ///////////////////
   // Check trigger
   ///////////////////
@@ -669,7 +673,7 @@ void AliTRDCalibTask::UserExec(Option_t *)
   Int_t nbtrackTPC = 0;
   
 
-  
+
   if (nbTracks <= 0.0) {
     
     if(fDebug > 1) {
@@ -680,8 +684,8 @@ void AliTRDCalibTask::UserExec(Option_t *)
     PostData(1, fListHist);
     return;
   }
-  
-  
+
+
   fESDfriend = dynamic_cast<AliESDfriend*> (fESD->FindListObject("AliESDfriend"));
   if(!fESDfriend){
     AliError("fESDfriend not available");
@@ -693,7 +697,7 @@ void AliTRDCalibTask::UserExec(Option_t *)
     PostData(1, fListHist);
     return;
   }
-  
+
   //printf("has friends\n");
 
   /////////////////////////////////////
@@ -759,7 +763,7 @@ void AliTRDCalibTask::UserExec(Option_t *)
       //printf("Not a good track\n");
       continue;
     }
-    
+
     // First Absolute gain calibration
     Int_t trdNTracklets = (Int_t) fkEsdTrack->GetTRDntracklets();
     Int_t trdNTrackletsPID = (Int_t) fkEsdTrack->GetTRDntrackletsPID(); 
@@ -784,9 +788,14 @@ void AliTRDCalibTask::UserExec(Option_t *)
       fTRDCalibraFillHisto->UpdateHistogramsV1(fTrdTrack,fkEsdTrack);
       //printf("Fill fTRDCalibraFillHisto\n");
     }
-      
-      
-         
+
+    const Double_t mag = AliTRDdEdxUtils::IsExBOn() ? fESD->GetMagneticField() : -1;
+    const Int_t charge = AliTRDdEdxUtils::IsExBOn() ? fkEsdTrack->Charge() : -1;
+    const Double_t toTPCscale = AliTRDdEdxUtils::GetCalibTPCscale(fkEsdTrack->GetTPCncls(), fkEsdTrack->GetTPCsignal());
+    if(toTPCscale>0){
+      AliTRDdEdxUtils::FillCalibHist(fTrdTrack, 0, mag, charge, toTPCscale);
+    }
+
     //////////////////////////////////
     // Debug 
     ////////////////////////////////
index c6ff2d3..677154e 100644 (file)
@@ -65,6 +65,7 @@
 #include "AliTRDCommonParam.h"
 #include "AliCDBManager.h"
 #include "AliCDBEntry.h"
+#include "AliTRDdEdxUtils.h"
 
 
 ClassImp(AliTRDPreprocessorOffline)
@@ -146,6 +147,8 @@ AliTRDPreprocessorOffline::~AliTRDPreprocessorOffline() {
   if(fPH2d) delete fPH2d;
   if(fPRF2d) delete fPRF2d;
   if(fSparse) delete fSparse;
+  AliTRDdEdxUtils::DeleteCalibHist();
+  AliTRDdEdxUtils::DeleteCalibObj();
   if(fAliTRDCalibraVdriftLinearFit) delete fAliTRDCalibraVdriftLinearFit;
   if(fAliTRDCalibraExbAltFit) delete fAliTRDCalibraExbAltFit;
   if(fNEvents) delete fNEvents;
@@ -167,8 +170,10 @@ void AliTRDPreprocessorOffline::Process(const Char_t* file, Int_t startRunNumber
     CalibGain(file,startRunNumber,endRunNumber,ocdbStorage);
     CalibChamberStatus(file,startRunNumber,endRunNumber,ocdbStorage);
     CalibExbAlt(file,startRunNumber,endRunNumber,ocdbStorage);
-    
+
   }
+
+  CalibPHQ(file, startRunNumber, endRunNumber, ocdbStorage);
   
   PrintStatus();
   
@@ -343,6 +348,27 @@ void AliTRDPreprocessorOffline::CalibPRF(const Char_t* file, Int_t startRunNumbe
   
 }
 //________________________________________________________________________________________________________________
+void AliTRDPreprocessorOffline::CalibPHQ(const Char_t* file, Int_t startRunNumber, Int_t endRunNumber, TString ocdbStorage)
+{
+  //
+  // make calibration of puls height Q
+  // Input parameters:
+  //      startRunNumber, endRunNumber     - run validity period 
+  //      ocdbStorage                      - path to the OCDB storage
+  //                                       - if empty - local storage 'pwd' uesed
+  //
+
+  if (ocdbStorage.Length()<=0) ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
+  //printf("test %s\n", ocdbStorage.Data());
+
+  if(!ReadPHQGlobal(file)) return;
+
+  if(!AnalyzePHQ(startRunNumber)) return;
+
+  UpdateOCDBPHQ(startRunNumber,endRunNumber,ocdbStorage);
+}
+
+//________________________________________________________________________________________________________________
 
 void AliTRDPreprocessorOffline::CalibChamberStatus(const Char_t* file, Int_t startRunNumber, Int_t endRunNumber, TString ocdbStorage){
   //
@@ -424,7 +450,7 @@ Bool_t AliTRDPreprocessorOffline::Init(const Char_t* fileName){
    
   if(fVersionVdriftUsed == 0) fStatusPos = fStatusPos |kVdriftErrorOld;
   if(fVersionGainUsed == 0) fStatusPos = fStatusPos | kGainErrorOld;
+
   return kTRUE;
   
 }
@@ -449,6 +475,17 @@ Bool_t AliTRDPreprocessorOffline::ReadStatusGlobal(const Char_t* fileName){
 }
 //___________________________________________________________________________________________________________________
 
+Bool_t AliTRDPreprocessorOffline::ReadPHQGlobal(const Char_t* fileName)
+{
+  //
+  // read calibration entries from file
+  //
+
+  return AliTRDdEdxUtils::ReadCalibHist(fileName, fNameList);
+}
+
+//___________________________________________________________________________________________________________________
+
 Bool_t AliTRDPreprocessorOffline::ReadGainGlobal(const Char_t* fileName){
   //
   // read calibration entries from file
@@ -894,6 +931,25 @@ Bool_t AliTRDPreprocessorOffline::AnalyzePRF(){
 }
 
 //_____________________________________________________________________________
+Bool_t AliTRDPreprocessorOffline::AnalyzePHQ(Int_t startRunNumber)
+{
+  //
+  //Produce PHQ calibration results
+  //
+  AliTRDdEdxUtils::PrintControl();
+
+  for(Int_t iter=0; iter<8; iter++){
+    THnSparseD *hi = (THnSparseD*) AliTRDdEdxUtils::GetHistPHQ()->At(iter);
+    TObjArray *obji = AliTRDdEdxUtils::GetCalibObj(hi, startRunNumber);
+    //printf("test analyze %s\n", obji->GetName());
+    AliTRDdEdxUtils::GetObjPHQ()->AddAt(obji, iter);
+  }
+
+  fCalibObjects->AddAt(AliTRDdEdxUtils::GetObjPHQ(), kPHQ);
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
 Bool_t AliTRDPreprocessorOffline::AnalyzeChamberStatus()
 {
   //
@@ -1255,6 +1311,27 @@ void AliTRDPreprocessorOffline::UpdateOCDBExBAlt(Int_t startRunNumber, Int_t end
 
 
  }
+//_________________________________________________________________________________________________________________
+void AliTRDPreprocessorOffline::UpdateOCDBPHQ(Int_t startRunNumber, Int_t endRunNumber, const Char_t *storagePath)
+{
+  //
+  // Update OCDB entry
+  //
+  AliCDBMetaData *metaData= new AliCDBMetaData();
+  metaData->SetObjectClassName("TObjArray");
+  metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
+  metaData->AddDateToComment();
+  metaData->SetBeamPeriod(1);
+
+  AliCDBId id1("TRD/Calib/PHQ", startRunNumber, endRunNumber);
+  AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
+  TObjArray *cobj = (TObjArray *) fCalibObjects->At(kPHQ);
+  if(cobj){
+    //cobj->Print();
+    gStorage->Put(cobj, id1, metaData);
+  }
+}
+
  //_________________________________________________________________________________________________________________
  void AliTRDPreprocessorOffline::UpdateOCDBChamberStatus(Int_t startRunNumber, Int_t endRunNumber, const Char_t *storagePath){
    //
index bfa3fb8..eab77e5 100644 (file)
@@ -36,7 +36,8 @@ public:
        kChamberStatus = 7,
        kPRF = 8,
        kExbAlt = 9,
-       kNumCalibObjs = 10
+        kPHQ = 10,
+       kNumCalibObjs = 11
   };
   enum { kGainNotEnoughStatsButFill = 2,
         kVdriftNotEnoughStatsButFill = 4,
@@ -151,6 +152,7 @@ public:
   void CalibGain(const Char_t* file, Int_t startRunNumber, Int_t endRunNumber,  TString  ocdbStorage="");
   void CalibPRF(const Char_t* file, Int_t startRunNumber, Int_t endRunNumber,  TString  ocdbStorage="");
   void CalibChamberStatus(const Char_t* file, Int_t startRunNumber, Int_t endRunNumber, TString ocdbStorage="");
+  void CalibPHQ(const Char_t* file, Int_t startRunNumber, Int_t endRunNumber, TString ocdbStorage);
 
   Bool_t ReadStatusGlobal(const Char_t* fileName="CalibObjects.root");
   Bool_t ReadGainGlobal(const Char_t* fileName="CalibObjects.root");
@@ -158,6 +160,7 @@ public:
   Bool_t ReadVdriftLinearFitGlobal(const Char_t* fileName="CalibObjects.root");
   Bool_t ReadExbAltFitGlobal(const Char_t* fileName="CalibObjects.root");
   Bool_t ReadPRFGlobal(const Char_t* fileName="CalibObjects.root");
+  Bool_t ReadPHQGlobal(const Char_t* fileName);
 
   Bool_t AnalyzeGain(); 
   Bool_t AnalyzeVdriftT0(); 
@@ -165,7 +168,8 @@ public:
   Bool_t AnalyzeExbAltFit();
   Bool_t AnalyzePRF();
   Bool_t AnalyzeChamberStatus(); 
-  
+  Bool_t AnalyzePHQ(Int_t startRunNumber);
+
   void CorrectFromDetGainUsed();
   void CorrectFromDetVdriftUsed();
   
@@ -176,6 +180,7 @@ public:
   void UpdateOCDBGain(Int_t  startRunNumber, Int_t endRunNumber, const char* storagePath);
   void UpdateOCDBPRF(Int_t  startRunNumber, Int_t endRunNumber, const char* storagePath);
   void UpdateOCDBChamberStatus(Int_t startRunNumber, Int_t endRunNumber, const Char_t *storagePath);
+  void UpdateOCDBPHQ(Int_t startRunNumber, Int_t endRunNumber, const Char_t *storagePath);
 
   Bool_t ValidateGain();
   Bool_t ValidateVdrift();
index 16fc2ac..1b229df 100644 (file)
@@ -321,6 +321,9 @@ const TObject *AliTRDcalibDB::GetCachedCDBObject(Int_t id)
     case kIDAttach : 
       return CacheCDBEntry(kIDAttach            ,"TRD/Calib/TrkAttach"); 
       break;
+    case kIDPHQ :
+      return CacheCDBEntry(kIDPHQ               ,"TRD/Calib/PHQ");
+      break;
   }
 
   return 0;
@@ -548,6 +551,16 @@ const AliTRDCalDet *AliTRDcalibDB::GetVdriftDet()
 }
 
 //_____________________________________________________________________________
+TObjArray * AliTRDcalibDB::GetPHQ()
+{
+  //
+  //return PHQ calibration object
+  //
+  TObjArray *arr = (TObjArray *) (GetCachedCDBObject(kIDPHQ));
+  return arr;
+}
+
+//_____________________________________________________________________________
 Float_t AliTRDcalibDB::GetVdriftAverage(Int_t det)
 {
   //
index 64619f1..debe520 100644 (file)
@@ -59,7 +59,7 @@ class AliTRDcalibDB : public TObject {
   Float_t                             GetVdriftAverage(Int_t det);
         AliTRDCalROC                 *GetVdriftROC(Int_t det);
   const AliTRDCalDet                 *GetVdriftDet();
-
+  TObjArray                          *GetPHQ();
   const AliTRDCalDet                 *GetExBDet();
 
   Float_t                             GetT0(Int_t det, Int_t col, Int_t row);
@@ -146,6 +146,7 @@ class AliTRDcalibDB : public TObject {
        , kIDPadStatus
        , kIDDCS
        , kIDAttach
+       , kIDPHQ
        , kCDBCacheSize };         // IDs of cached objects
 
   const TObject *GetCachedCDBObject(Int_t id);
index a19939a..2c7981f 100644 (file)
@@ -16,6 +16,7 @@
 
 /* $Id$ */
 
+#include "TVectorT.h"
 #include "AliLog.h"
 #include "AliESDtrack.h"
 #include "AliTracker.h"
@@ -26,6 +27,7 @@
 #include "AliTRDReconstructor.h"
 #include "AliTRDPIDResponse.h"
 #include "AliTRDrecoParam.h"
+#include "AliTRDdEdxUtils.h"
 
 ClassImp(AliTRDtrackV1)
 
@@ -45,6 +47,7 @@ AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
   ,fStatus(0)
   ,fESDid(0)
   ,fDE(0.)
+  ,fTruncatedMean(0)
   ,fkReconstructor(NULL)
   ,fBackupTrack(NULL)
   ,fTrackLow(NULL)
@@ -71,6 +74,7 @@ AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
   ,fStatus(ref.fStatus)
   ,fESDid(ref.fESDid)
   ,fDE(ref.fDE)
+  ,fTruncatedMean(ref.fTruncatedMean)
   ,fkReconstructor(ref.fkReconstructor)
   ,fBackupTrack(NULL)
   ,fTrackLow(NULL)
@@ -101,6 +105,7 @@ AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
   ,fStatus(0)
   ,fESDid(0)
   ,fDE(0.)
+  ,fTruncatedMean(0)
   ,fkReconstructor(NULL)
   ,fBackupTrack(NULL)
   ,fTrackLow(NULL)
@@ -154,6 +159,7 @@ AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 * const trklts, const Double_t p[5], c
   ,fStatus(0)
   ,fESDid(0)
   ,fDE(0.)
+  ,fTruncatedMean(0)
   ,fkReconstructor(NULL)
   ,fBackupTrack(NULL)
   ,fTrackLow(NULL)
@@ -381,6 +387,14 @@ Bool_t AliTRDtrackV1::CookPID()
     }
   }
   pidResponse->GetResponse(nslices, dEdx, trackletP, fPID);
+
+  //do truncated mean
+  //ncls needs to be included!! todo!!
+  AliTRDdEdxUtils::SetObjPHQ(AliTRDcalibDB::Instance()->GetPHQ());
+  const Double_t mag    = AliTRDdEdxUtils::IsExBOn() ? GetBz()  : -1;
+  const Double_t charge = AliTRDdEdxUtils::IsExBOn() ? Charge() : -1;
+  fTruncatedMean = CookTruncatedMean(0, mag, charge, kTRUE);
+
   return kTRUE;
 }
 
@@ -907,4 +921,32 @@ void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
   }
   // store PID probabilities
   track->SetTRDpid(fPID);
+
+  //store truncated mean
+  track->SetTRDsignal(fTruncatedMean);
 }
+
+//_______________________________________________________________
+Double_t  AliTRDtrackV1::CookTruncatedMean(const Bool_t kinvq, const Double_t mag, const Int_t charge, const Int_t kcalib, TVectorD *Qs, TVectorD *Xs, Int_t timeBin0, Int_t timeBin1, Int_t tstep) const
+{
+  //
+  //Origin: Xianguo Lu <xianguo.lu@cern.ch>, Marian Ivanov <marian.ivanov@cern.ch>
+  //
+
+  TVectorD arrayQ(200), arrayX(200);
+  Int_t ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, this, timeBin0, timeBin1, tstep);
+
+  const TObjArray *cobj = kcalib ? AliTRDdEdxUtils::GetObjPHQ(kinvq, mag, charge) : NULL;
+  
+  const Double_t tmean = AliTRDdEdxUtils::ToyCook(kinvq, ncls, &arrayQ, &arrayX, cobj);
+
+  const Int_t nch = AliTRDdEdxUtils::UpdateArrayX(ncls, &arrayX);
+
+  if(Qs && Xs){
+    (*Qs)=arrayQ;
+    (*Xs)=arrayX;
+  }
+
+  return AliTRDdEdxUtils::GetSignal(nch, ncls, tmean);
+}
+
index 4f2b154..630423f 100644 (file)
@@ -19,6 +19,8 @@
 #include "AliTRDseedV1.h"
 #endif
 
+template <typename Value> class TVectorT;
+typedef class TVectorT<Double_t> TVectorD;
 class AliESDtrack;
 class AliTRDcluster;
 class AliTRDReconstructor;
@@ -76,6 +78,8 @@ public:
   virtual void   Copy(TObject &ref) const;
  
   Bool_t         CookPID();
+  Double_t       CookTruncatedMean(const Bool_t kinvq, const Double_t mag, const Int_t charge, const Int_t kcalib, TVectorD *Qs=NULL, TVectorD *Xs=NULL, Int_t timeBin0=-1, Int_t timeBin1=1000, Int_t tstep=1) const;
   Int_t          CookLabel(Float_t wrong, Int_t *labs=NULL, Float_t *freq=NULL);
   AliTRDtrackV1* GetBackupTrack() const {return fBackupTrack;}
   Double_t       GetBudget(Int_t i) const { return fBudget[i];}
@@ -145,6 +149,8 @@ private:
   Double32_t   fPID[AliPID::kSPECIES]; //  PID probabilities
   Double32_t   fBudget[3];             //  Integrated material budget
   Double32_t   fDE;                    //  Integrated delta energy
+  Double32_t   fTruncatedMean;         //  Truncated mean
+
   const AliTRDReconstructor *fkReconstructor;//! reconstructor link 
   AliTRDtrackV1 *fBackupTrack;         //! Backup track
   AliTRDseedV1  *fTracklet[kNplane];   //  Tracklets array defining the track
index d2514fa..6ca5528 100644 (file)
@@ -33,6 +33,7 @@ set ( SRCS
     AliTRDtransform.cxx
     AliTRDtrackletOflHelper.cxx
     AliTRDpidUtil.cxx
+    AliTRDdEdxUtils.cxx
     AliTRDpidESD.cxx
     AliTRDReconstructor.cxx
     AliTRDseedV1.cxx AliTRDtrackV1.cxx
index 9379eb0..be3c58e 100644 (file)
@@ -14,6 +14,7 @@
 
 #pragma link C++ class  AliTRDpidESD+;
 #pragma link C++ class  AliTRDpidUtil+;
+#pragma link C++ class  AliTRDdEdxUtils+;
 
 #pragma link C++ class  AliTRDReconstructor+;