From e2a1c98b27a8338047c9de092a4cb0c977e50954 Mon Sep 17 00:00:00 2001 From: abercuci Date: Mon, 26 Mar 2012 11:21:32 +0000 Subject: [PATCH] Truncated mean PID method by Xianguo --- TRD/AliTRDCalibTask.cxx | 27 +++++++---- TRD/AliTRDPreprocessorOffline.cxx | 81 ++++++++++++++++++++++++++++++- TRD/AliTRDPreprocessorOffline.h | 9 +++- TRD/AliTRDcalibDB.cxx | 13 +++++ TRD/AliTRDcalibDB.h | 3 +- TRD/AliTRDtrackV1.cxx | 42 ++++++++++++++++ TRD/AliTRDtrackV1.h | 6 +++ TRD/CMakelibTRDrec.pkg | 1 + TRD/TRDrecLinkDef.h | 1 + 9 files changed, 169 insertions(+), 14 deletions(-) diff --git a/TRD/AliTRDCalibTask.cxx b/TRD/AliTRDCalibTask.cxx index 52aa79fb398..6df00b04cf0 100644 --- a/TRD/AliTRDCalibTask.cxx +++ b/TRD/AliTRDCalibTask.cxx @@ -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 (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 //////////////////////////////// diff --git a/TRD/AliTRDPreprocessorOffline.cxx b/TRD/AliTRDPreprocessorOffline.cxx index c6ff2d393f8..677154e697b 100644 --- a/TRD/AliTRDPreprocessorOffline.cxx +++ b/TRD/AliTRDPreprocessorOffline.cxx @@ -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(); @@ -342,6 +347,27 @@ void AliTRDPreprocessorOffline::CalibPRF(const Char_t* file, Int_t startRunNumbe UpdateOCDBPRF(startRunNumber,endRunNumber,ocdbStorage); } +//________________________________________________________________________________________________________________ +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 @@ -893,6 +930,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){ // diff --git a/TRD/AliTRDPreprocessorOffline.h b/TRD/AliTRDPreprocessorOffline.h index bfa3fb812c3..eab77e5a1f1 100644 --- a/TRD/AliTRDPreprocessorOffline.h +++ b/TRD/AliTRDPreprocessorOffline.h @@ -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(); diff --git a/TRD/AliTRDcalibDB.cxx b/TRD/AliTRDcalibDB.cxx index 16fc2acde73..1b229dfff96 100644 --- a/TRD/AliTRDcalibDB.cxx +++ b/TRD/AliTRDcalibDB.cxx @@ -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; @@ -547,6 +550,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) { diff --git a/TRD/AliTRDcalibDB.h b/TRD/AliTRDcalibDB.h index 64619f19dc0..debe5205b51 100644 --- a/TRD/AliTRDcalibDB.h +++ b/TRD/AliTRDcalibDB.h @@ -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); diff --git a/TRD/AliTRDtrackV1.cxx b/TRD/AliTRDtrackV1.cxx index a19939a2e29..2c7981fe80a 100644 --- a/TRD/AliTRDtrackV1.cxx +++ b/TRD/AliTRDtrackV1.cxx @@ -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 , Marian Ivanov + // + + 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); +} + diff --git a/TRD/AliTRDtrackV1.h b/TRD/AliTRDtrackV1.h index 4f2b1548071..630423fb69e 100644 --- a/TRD/AliTRDtrackV1.h +++ b/TRD/AliTRDtrackV1.h @@ -19,6 +19,8 @@ #include "AliTRDseedV1.h" #endif +template class TVectorT; +typedef class TVectorT 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 diff --git a/TRD/CMakelibTRDrec.pkg b/TRD/CMakelibTRDrec.pkg index d2514fab3a5..6ca55284ea7 100644 --- a/TRD/CMakelibTRDrec.pkg +++ b/TRD/CMakelibTRDrec.pkg @@ -33,6 +33,7 @@ set ( SRCS AliTRDtransform.cxx AliTRDtrackletOflHelper.cxx AliTRDpidUtil.cxx + AliTRDdEdxUtils.cxx AliTRDpidESD.cxx AliTRDReconstructor.cxx AliTRDseedV1.cxx AliTRDtrackV1.cxx diff --git a/TRD/TRDrecLinkDef.h b/TRD/TRDrecLinkDef.h index 9379eb07b25..be3c58e91be 100644 --- a/TRD/TRDrecLinkDef.h +++ b/TRD/TRDrecLinkDef.h @@ -14,6 +14,7 @@ #pragma link C++ class AliTRDpidESD+; #pragma link C++ class AliTRDpidUtil+; +#pragma link C++ class AliTRDdEdxUtils+; #pragma link C++ class AliTRDReconstructor+; -- 2.39.3