{
track -> SetReconstructor(fReconstructor);
-
- fReconstructor -> SetOption("nn");
+ (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork();
track -> CookPID();
if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion) + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
}
if(!IsInRange(momentum)) return 0x0;
- fReconstructor -> SetOption("!nn");
+ (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(/*kFALSE*/);
cTrack.CookPID();
if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
}
if(!IsInRange(momentum)) return 0x0;
- fReconstructor -> SetOption("nn");
+ (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork();
cTrack.CookPID();
if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
if(!IsInRange(momentum)) return 0x0;
- fReconstructor -> SetOption("!nn");
+ (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(/*kFALSE*/);
Float_t SumdEdx = 0;
Int_t iBin = FindBin(species, momentum);
AliTRDseedV1 *tracklet = 0x0;
}
if(!IsInRange(momentum)) return 0x0;
- fReconstructor -> SetOption("!nn");
+ (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(/*kFALSE*/);
Int_t iMomBin = fMomentumAxis->FindBin(momentum);
Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
Float_t *fdEdx;
ClassImp(AliTRDpidRefMaker)
//________________________________________________________________________
-AliTRDpidRefMaker::AliTRDpidRefMaker(UChar_t n, const char *name, const char *title)
+AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title)
:AliTRDrecoTask(name, title)
,fReconstructor(0x0)
,fV0s(0x0)
,fData(0x0)
,fRefPID(kMC)
,fRefP(kMC)
- ,fNslices(n)
,fTrainFreq(1.)
,fTestFreq(0.)
- ,fScale(1.)
,fLy(-1)
,fP(-1.)
{
OpenFile(0, "RECREATE");
fContainer = new TObjArray();
+ fContainer->SetName(Form("Moni%s", GetName()));
fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
-
- // open reference TTree
- OpenFile(1, "RECREATE");
- fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
- fData->Branch("l" , &fLy , "l/b");
- fData->Branch("p" , &fP , "p/F");
- fData->Branch("pid", fPID, Form("fPID[%d]/F", AliPID::kSPECIES));
- fData->Branch("dEdx" , fdEdx , Form("fdEdx[%d]/F", fNslices));
}
//________________________________________________________________________
}
// set reconstructor
- Float_t *dedx;
+ //Float_t *dedx;
TRDtrack -> SetReconstructor(fReconstructor);
// fill the momentum and dE/dx information
for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
- TRDtracklet->CookdEdx(fNslices);
- dedx = const_cast<Float_t*>(TRDtracklet->GetdEdx());
- for(Int_t is=fNslices; is--;) fdEdx[is] = dedx[is]/fScale;
-
+ if(!GetdEdx(TRDtracklet)) continue;
switch(fRefP){
case kMC:
if(!HasMCdata()){
return;
}
fLy = ily;
- fData->Fill();
+ Fill();
}
for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
}
+//________________________________________________________________________
+void AliTRDpidRefMaker::Fill()
+{
+ fData->Fill();
+}
//________________________________________________________________________
void AliTRDpidRefMaker::Terminate(Option_t *)
}
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::MakeTrainingList()
-{
- //
- // build the training lists for the neural networks
- //
-
- if (!fData) LoadFile(Form("TRD.%s.root", GetName()));
- if (!fData) {
- AliError("Tree for training list not available");
- return;
- }
-
- Int_t nPart[AliTRDCalPID::kNMom][AliTRDgeometry::kNlayer][AliPID::kSPECIES];
- memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*AliTRDgeometry::kNlayer*sizeof(Int_t));
-
- // set needed branches
- fData -> SetBranchAddress("l", &fLy);
- fData -> SetBranchAddress("dEdx", &fdEdx);
- fData -> SetBranchAddress("pid", &fPID);
- fData -> SetBranchAddress("p", &fP);
-
- AliTRDpidUtil *util = new AliTRDpidUtil();
-
- // start first loop to check total number of each particle type
- Int_t n(0);
- for(Int_t iEv=0; iEv < fData -> GetEntries(); iEv++){
- fData -> GetEntry(iEv);
-
- // get particle type
- Int_t pidPos = TMath::LocMax(AliPID::kSPECIES, fPID);
-
- // get momentum bin
- Int_t ip = util->GetMomentumBin(fP);
- if(ip<0){
- AliWarning(Form("Wrong momentum %f.", fP));
- continue;
- }
- nPart[ip][fLy][pidPos]++;
- n++;
- }
-
- AliDebug(2, Form("Particle multiplicities [%d]:", n));
- n=0;
- for(Int_t ip = 0; ip <AliTRDCalPID::kNMom; ip++){
- printf(" P[%2d] ", ip);
- for(Int_t is = 0; is <AliPID::kSPECIES; is++){
- printf("%s[", AliPID::ParticleShortName(is));
- for(Int_t il = 0; il <AliTRDgeometry::kNlayer; il++){
- printf("%d ", nPart[ip][il][is]);
- n+=nPart[ip][il][is];
- }
- printf("] ");
- }
- printf("\n");
- }
- AliDebug(2, Form("Particle multiplicities check [%d]:", n));
-
-
- // set training/test sample size per momentum interval
- Int_t iTrain[AliTRDCalPID::kNMom], iTest[AliTRDCalPID::kNMom];
- for(Int_t ip = 0; ip < AliTRDCalPID::kNMom; ip++){
- Int_t min = 1000000;//TMath::MinElement(AliPID::kSPECIES, nPart[ip][5]);
- AliDebug(10, Form("Ref Stat in pBin[%d]=%d", ip, min));
- iTrain[ip] = Int_t(min * fTrainFreq);
- iTest[ip] = Int_t(min * fTestFreq);
- AliDebug(10, Form("P_bin[%2d] Train[%d] Test[%d]", ip, iTrain[ip], iTest[ip]));
- }
-
-
- // start second loop to set the event lists
- // reset couters
- memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*AliTRDgeometry::kNlayer*sizeof(Int_t));
- for(Int_t iEv = 0; iEv < fData->GetEntries(); iEv++){
- fData->GetEntry(iEv);
- // get PID position
- Int_t pidPos = TMath::LocMax(AliPID::kSPECIES, fPID);
-
-
- Int_t iMomBin = util->GetMomentumBin(fP);
- if(nPart[iMomBin][fLy][pidPos] < iTrain[iMomBin]){
- fTrain[iMomBin][fLy]->Enter(iEv);
- nPart[iMomBin][fLy][pidPos]++;
- } else if(nPart[iMomBin][fLy][pidPos] < iTest[iMomBin]+iTrain[iMomBin]){
- fTest[iMomBin][fLy]->Enter(iEv);
- nPart[iMomBin][fLy][pidPos]++;
- } else continue;
- }
-
- AliDebug(2, "Particle multiplicities in both lists:");
- for(Int_t ip = 0; ip <AliTRDCalPID::kNMom; ip++){
- printf(" P[%2d] ", ip);
- for(Int_t is = 0; is <AliPID::kSPECIES; is++){
- Int_t m(0);
- for(Int_t il = 0; il <AliTRDgeometry::kNlayer; il++) m+=nPart[ip][il][is];
- printf("%s[%4d] ", AliPID::ParticleShortName(is), m);
- }
- printf("\n");
- }
-
- util->Delete();
- //delete util;
-}
-
-
-
-
//________________________________________________________________________
void AliTRDpidRefMaker::LoadFile(const Char_t *InFile)
{
class TObjArray;
class TEventList;
class AliTRDReconstructor;
+class AliTRDseedV1;
class AliTRDpidRefMaker : public AliTRDrecoTask
{
,kMC = 1 // use MC as reference
,kRec= 2 // use Reconstructed PID as reference
};
- AliTRDpidRefMaker(UChar_t n, const char *name=0, const char *title=0);
+ AliTRDpidRefMaker(const char *name=0, const char *title=0);
virtual ~AliTRDpidRefMaker();
void ConnectInputData(Option_t *opt);
void CreateOutputObjects();
void Exec(Option_t *option);
- Float_t GetScaledEdx() const { return fScale;}
+
void LoadContainer(const Char_t *InFileCont);
void LoadFile(const Char_t *InFile);
- void MakeTrainingList(); // build the training and the test list
void SetAbundance(Float_t train, Float_t test);
void SetRefPID(void *source, Float_t *pid);
void SetSource(ETRDpidRefMakerSource pid, ETRDpidRefMakerSource momentum) {fRefPID = pid; fRefP = momentum;}
- void SetScaledEdx(Float_t s) {fScale = s;}
+
void Terminate(Option_t *);
protected:
- virtual void MakeRefs(Int_t mombin)=0; // build the reference for a given momentum bin
+ virtual Float_t* GetdEdx(AliTRDseedV1*) = 0;
+ virtual Int_t GetNslices() = 0;
+ virtual void Fill();
+
AliTRDReconstructor *fReconstructor; //! reconstructor needed for recalculation the PID
TObjArray *fV0s; //! v0 array
TTree *fData; //! dEdx-P data
TEventList *fTest[AliTRDCalPID::kNMom][AliTRDgeometry::kNlayer]; // Test list for each momentum
ETRDpidRefMakerSource fRefPID; // reference PID source
ETRDpidRefMakerSource fRefP; // reference momentum source
- UChar_t fNslices; // number of dEdx slices used by implementation
Float_t fTrainFreq; //! training sample relative abundance
Float_t fTestFreq; //! testing sample relative abundance
- Float_t fScale; //! dEdx scaling
UChar_t fLy; //! TRD layer
Float_t fP; //! momentum
Float_t fdEdx[10]; //! dEdx array
//
///////////////////////////////////////////////////////////////////////////////
-#include <TROOT.h>
#include <TFile.h>
+#include <TROOT.h>
#include <TTree.h>
+#include <TMath.h>
+#include <TEventList.h>
#include <TH2D.h>
#include <TH2I.h>
-#include <TParticle.h>
-#include <TParticle.h>
#include <TPrincipal.h>
#include <TVector3.h>
#include <TLinearFitter.h>
#include <TEllipse.h>
#include <TMarker.h>
-#include "AliPID.h"
-#include "AliESD.h"
-#include "AliRun.h"
-#include "AliRunLoader.h"
-#include "AliStack.h"
-#include "AliHeader.h"
-#include "AliGenEventHeader.h"
-#include "AliESDtrack.h"
-
+#include "AliLog.h"
+#include "../../STAT/TKDPDF.h"
#include "AliTRDpidRefMakerLQ.h"
#include "../Cal/AliTRDCalPID.h"
+#include "AliTRDseedV1.h"
#include "AliTRDcalibDB.h"
#include "AliTRDgeometry.h"
ClassImp(AliTRDpidRefMakerLQ)
-TLinearFitter *AliTRDpidRefMakerLQ::fgFitter2D2 = 0x0;
-TLinearFitter *AliTRDpidRefMakerLQ::fgFitter2D1 = 0x0;
-
//__________________________________________________________________
AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
- :AliTRDpidRefMaker(UChar_t(AliTRDpidUtil::kLQslices), "PidRefMakerLQ", "PID(LQ) Reference Maker")
+ :AliTRDpidRefMaker("PidRefMakerLQ", "PID(LQ) Reference Maker")
+ ,fPbin(-1)
+ ,fSbin(-1)
{
//
// AliTRDpidRefMakerLQ default constructor
//
- // histogram settings
- for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
- fH2dEdx[ispec] = 0x0;
- fPrinc[ispec] = 0x0;
- }
+ memset(fH2dEdx, 0x0, AliPID::kSPECIES*sizeof(TH2*));
}
//__________________________________________________________________
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
if(fH2dEdx[ispec]) delete fH2dEdx[ispec];
- if(fPrinc[ispec]) delete fPrinc[ispec];
}
- if(fgFitter2D1){ delete fgFitter2D1; fgFitter2D1 = 0x0;}
- if(fgFitter2D2){ delete fgFitter2D2; fgFitter2D2 = 0x0;}
-
}
-
-//__________________________________________________________________
-void AliTRDpidRefMakerLQ::MakeRefs(Int_t /*pbin*/)
+//________________________________________________________________________
+void AliTRDpidRefMakerLQ::CreateOutputObjects()
{
- // Build, Fill and write to file the histograms used for PID.
- // The simulations are looked in the
- // directories with the general form Form("p%3.1f", momentum)
- // starting from dir (default .). Here momentum belongs to the list
- // of known momentum to PID (see trackMomentum).
- // The output histograms are
- // written to the file "File" in cwd (default
- // TRDPIDHistograms.root). In order to build a DB entry
- // consider running $ALICE_ROOT/Cal/AliTRDpidCDB.C
- //
- // Author:
- // Alex Bercuci (A.Bercuci@gsi.de)
-
- Int_t partCode[AliPID::kSPECIES] =
- {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
-
- // check and retrive number of directories in the production
- const char *dir=".", *File="LQ.root";
- Int_t nBatches;
- if(!(nBatches = CheckProdDirTree(dir))) return;
-
-
- // Number of Time bins
- Int_t nTimeBins;
- AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
- if(!calibration){
- AliError("No AliTRDcalibDB available.");
- return;
- } else {
- nTimeBins = calibration->GetNumberOfTimeBins();
-// if (calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
-// else {
-// AliError("No run number set.");
-// return kFALSE;
-// }
- }
-
- // Build PID reference/working objects
- fH1TB[0] = new TH1F("h1TB_0", "", nTimeBins, -.5, nTimeBins-.5);
- fH1TB[0]->SetLineColor(4);
- fH1TB[1] = new TH1F("h1TB_1", "", nTimeBins, -.5, nTimeBins-.5);
- fH1TB[1]->SetLineColor(2);
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) if(!fPrinc[ispec]) fPrinc[ispec] = new TPrincipal(2, "ND");
-
- // Print statistics header
- Int_t nPart[AliPID::kSPECIES], nTotPart;
- printf("P[GeV/c] ");
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPID::GetPartSymb(ispec));
- printf("\n-----------------------------------------------\n");
-
- Float_t trackMomentum[AliTRDCalPID::kNMom];
- //Float_t trackSegLength[AliTRDCalPID::kNLength];
- for(int i=0; i<AliTRDCalPID::kNMom; i++) trackMomentum[i] = AliTRDCalPID::GetMomentum(i);
- AliRunLoader *fRunLoader = 0x0;
- TFile *esdFile = 0x0;
- TTree *esdTree = 0x0;
- AliESD *esd = 0x0;
- AliESDtrack *esdTrack = 0x0;
-
- //
- // Momentum loop
- for (Int_t imom = 0; imom < AliTRDCalPID::kNMom; imom++) {
- Reset();
-
- // init statistics for momentum
- nTotPart = 0;
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
-
- // loop over production directories
- for(Int_t ibatch = 0; ibatch<nBatches; ibatch++){
- // open run loader and load gAlice, kinematics and header
- fRunLoader = AliRunLoader::Open(Form("%s/%3.1fGeV/%03d/galice.root", dir, trackMomentum[imom], ibatch));
- if (!fRunLoader) {
- AliError(Form("Getting run loader for momentum %3.1f GeV/c batch %03d failed.", trackMomentum[imom], ibatch));
- return;
- }
- TString s; s.Form("%s/%3.1fGeV/%03d/", dir, trackMomentum[imom], ibatch);
- fRunLoader->SetDirName(s);
- fRunLoader->LoadgAlice();
- gAlice = fRunLoader->GetAliRun();
- if (!gAlice) {
- AliError(Form("galice object not found for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch));
- return;
- }
- fRunLoader->LoadKinematics();
- fRunLoader->LoadHeader();
-
- // open the ESD file
- esdFile = TFile::Open(Form("%s/%3.1fGeV/%03d/AliESDs.root", dir, trackMomentum[imom], ibatch));
- if (!esdFile || esdFile->IsZombie()) {
- AliError(Form("Opening ESD file failed for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch));
- return;
- }
- esdTree = (TTree*)esdFile->Get("esdTree");
- if (!esdTree) {
- AliError(Form("ESD tree not found for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch));
- return;
- }
- esd = new AliESD;
- esdTree->SetBranchAddress("ESD", &esd);
-
-
- // Event loop
- for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
-
- // Load MC info
- fRunLoader->GetEvent(iEvent);
- AliStack* stack = AliRunLoader::Instance()->Stack();
- TArrayF vertex(3);
- fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
-
- // Load event summary data
- esdTree->GetEvent(iEvent);
- if (!esd) {
- AliWarning(Form("ESD object not found for event %d. (@ momentum %3.1f GeV/c, batch %03d)", iEvent, trackMomentum[imom], ibatch));
- continue;
- }
-
- // Track loop
- for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
- esdTrack = esd->GetTrack(iTrack);
-
- //if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
-
- if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
- if(esdTrack->GetConstrainedChi2() > 1E9) continue;
- if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
- if (esdTrack->GetTRDsignal() == 0.) continue;
-
- // read MC info
- Int_t label = esdTrack->GetLabel();
- if(label<0) continue;
- if (label > stack->GetNtrack()) continue; // background
- TParticle* particle = stack->Particle(label);
- if(!particle){
- AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ momentum %3.1f batch %03d event %d track %d]", label, trackMomentum[imom], ibatch, iEvent, iTrack));
- continue;
- }
- if(particle->Pt() < 1.E-3) continue;
- // if (TMath::Abs(particle->Eta()) > 0.3) continue;
- TVector3 dVertex(particle->Vx() - vertex[0],
- particle->Vy() - vertex[1],
- particle->Vz() - vertex[2]);
- if (dVertex.Mag() > 1.E-4){
- //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
- //particle->Print();
- continue;
- }
- Int_t iGen = -1;
- for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
- if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
- iGen = ispec;
- break;
- }
- if(iGen<0) continue;
-
- nPart[iGen]++; nTotPart++;
-
- Float_t mom;
- //Float_t length;
- Double_t dedx[AliTRDCalPID::kNSlicesLQ], dEdx;
- Int_t timebin;
- for (Int_t iLayer=0; iLayer<AliTRDgeometry::kNlayer; iLayer++){
- // read data for track segment
- for(int iSlice=0; iSlice<AliTRDCalPID::kNSlicesLQ; iSlice++)
- dedx[iSlice] = esdTrack->GetTRDslice(iLayer, iSlice);
- dEdx = esdTrack->GetTRDslice(iLayer, -1);
- timebin = esdTrack->GetTRDTimBin(iLayer);
-
- // check data
- if ((dEdx <= 0.) || (timebin <= -1.)) continue;
-
- // retrive kinematic info for this track segment
- //if(!AliTRDpidESD::RecalculateTrackSegmentKine(esdTrack, iLayer, mom, length)) continue;
- mom = esdTrack->GetOuterParam()->GetP();
-
- // find segment length and momentum bin
- Int_t jmom = 1, refMom = -1;
- while(jmom<AliTRDCalPID::kNMom-1 && mom>trackMomentum[jmom]) jmom++;
- if(TMath::Abs(trackMomentum[jmom-1] - mom) < trackMomentum[jmom-1] * .2) refMom = jmom-1;
- else if(TMath::Abs(trackMomentum[jmom] - mom) < trackMomentum[jmom] * .2) refMom = jmom;
- if(refMom<0){
- AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ momentum %3.1f batch %03d event %d track %d]", iLayer, trackMomentum[imom], ibatch, iEvent, iTrack));
- continue;
- }
- /*while(jleng<AliTRDCalPID::kNLength-1 && length>trackSegLength[jleng]) jleng++;*/
-
- // this track segment has fulfilled all requierments
- //nPlanePID++;
-
- if(dedx[0] > 0. && dedx[1] > 0.){
- dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
- fPrinc[iGen]->AddRow(dedx);
- }
- fH1TB[(iGen>0)?1:0]->Fill(timebin);
- } // end plane loop
- } // end track loop
- } // end events loop
-
- delete esd; esd = 0x0;
- esdFile->Close();
- delete esdFile; esdFile = 0x0;
-
- fRunLoader->UnloadHeader();
- fRunLoader->UnloadKinematics();
- delete fRunLoader; fRunLoader = 0x0;
- } // end directory loop
-
- // use data to prepare references
- Prepare2D();
- // save references
- SaveReferences(imom, File);
-
-
- // print momentum statistics
- printf(" %3.1f ", trackMomentum[imom]);
- for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
- printf("\n");
- } // end momentum loop
-
- TFile *fSave = 0x0;
- TListIter it((TList*)gROOT->GetListOfFiles());
- while((fSave=(TFile*)it.Next()))
- if(strcmp(File, fSave->GetName())==0) break;
+ // Create histograms
+ // Called once
- fSave->cd();
- fSave->Close();
- delete fSave;
+ AliTRDpidRefMaker::CreateOutputObjects();
- return;
+ fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
+ fData->Branch("s", &fSbin, "l/b");
+ fData->Branch("p", &fPbin, "p/b");
+ fData->Branch("dEdx" , fdEdx , Form("dEdx[%d]/F", GetNslices()));
}
-//__________________________________________________________________
-Double_t AliTRDpidRefMakerLQ::Estimate2D2(TH2 * const h, Float_t &x, Float_t &y)
-{
- //
- // Linear interpolation of data point with a parabolic expresion using
- // the logarithm of the bin content in a close neighbourhood. It is
- // assumed that the bin content of h is in number of events !
- //
- // Observation:
- // This function have to be used when the statistics of each bin is
- // sufficient and all bins are populated. For cases of sparse data
- // please refere to Estimate2D1().
- //
- // Author : Alex Bercuci (A.Bercuci@gsi.de)
- //
-
- if(!h){
- AliErrorGeneral("AliTRDpidRefMakerLQ::Estimate2D2()", "No histogram defined.");
- return 0.;
- }
-
- TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
- Int_t binx = ax->FindBin(x);
- Int_t biny = ay->FindBin(y);
- Int_t nbinsx = ax->GetNbins();
- Int_t nbinsy = ay->GetNbins();
- Double_t p[2];
- Double_t entries;
-
- // late construction of fitter
- if(!fgFitter2D2) fgFitter2D2 = new TLinearFitter(6, "1++x++y++x*x++y*y++x*y");
-
- fgFitter2D2->ClearPoints();
- Int_t npoints=0;
- Int_t binx0, binx1, biny0, biny1;
- for(int bin=0; bin<5; bin++){
- binx0 = TMath::Max(1, binx-bin);
- binx1 = TMath::Min(nbinsx, binx+bin);
- for(int ibin=binx0; ibin<=binx1; ibin++){
- biny0 = TMath::Max(1, biny-bin);
- biny1 = TMath::Min(nbinsy, biny+bin);
- for(int jbin=biny0; jbin<=biny1; jbin++){
- if(ibin != binx0 && ibin != binx1 && jbin != biny0 && jbin != biny1) continue;
- if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
- p[0] = ax->GetBinCenter(ibin);
- p[1] = ay->GetBinCenter(jbin);
- fgFitter2D2->AddPoint(p, log(entries), 1./sqrt(entries));
- npoints++;
- }
- }
- if(npoints>=25) break;
- }
- if(fgFitter2D2->Eval() == 1){
- printf("<I2> x = %9.4f y = %9.4f\n", x, y);
- printf("\tbinx %d biny %d\n", binx, biny);
- printf("\tpoints %d\n", npoints);
-
- return 0.;
- }
- TVectorD vec(6);
- fgFitter2D2->GetParameters(vec);
- Double_t result = vec[0] + x*vec[1] + y*vec[2] + x*x*vec[3] + y*y*vec[4] + x*y*vec[5];
- return exp(result);
-}
-
-//__________________________________________________________________
-Double_t AliTRDpidRefMakerLQ::Estimate2D1(TH2 * const h, Float_t &x, Float_t &y
- , const Float_t &dCT
- , const Float_t &rmin
- , const Float_t &rmax)
+//________________________________________________________________________
+Float_t* AliTRDpidRefMakerLQ::GetdEdx(AliTRDseedV1 *trklt)
{
- //
- // Linear interpolation of data point with a plane using
- // the logarithm of the bin content in the area defined by the
- // d(cos(phi)) and dr=(rmin, rmax). It is assumed that the bin content
- // of h is number of events !
- //
-
- if(!h){
- AliErrorGeneral("AliTRDpidRefMakerLQ::Estimate2D1()", "No histogram defined.");
- return 0.;
- }
-
- TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
-// Int_t binx = ax->FindBin(x);
-// Int_t biny = ay->FindBin(y);
- Int_t nbinsx = ax->GetNbins();
- Int_t nbinsy = ay->GetNbins();
- Double_t p[2];
- Double_t entries;
- Double_t rxy = sqrt(x*x + y*y), rpxy;
-
- // late construction of fitter
- if(!fgFitter2D1) fgFitter2D1 = new TLinearFitter(3, "1++x++y");
-
- fgFitter2D1->ClearPoints();
- Int_t npoints=0;
- for(int ibin=1; ibin<=nbinsx; ibin++){
- for(int jbin=1; jbin<=nbinsy; jbin++){
- if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
- p[0] = ax->GetBinCenter(ibin);
- p[1] = ay->GetBinCenter(jbin);
- rpxy = sqrt(p[0]*p[0] + p[1]*p[1]);
- if((x*p[0] + y*p[1])/rxy/rpxy < dCT) continue;
- if(rpxy<rmin || rpxy > rmax) continue;
-
- fgFitter2D1->AddPoint(p, log(entries), 1./sqrt(entries));
- npoints++;
- }
- }
- if(npoints<15) return 0.;
- if(fgFitter2D1->Eval() == 1){
- printf("<O2> x = %9.4f y = %9.4f\n", x, y);
- printf("\tpoints %d\n", npoints);
- return 0.;
- }
- TVectorD vec(3);
- fgFitter2D1->GetParameters(vec);
- Double_t result = vec[0] + x*vec[1] + y*vec[2];
- return exp(result);
+ trklt->CookdEdx(AliTRDpidUtil::kLQslices);
+ const Float_t *dedx = trklt->GetdEdx();
+ if(dedx[0]+dedx[1] <= 0.) return 0x0;
+ if(dedx[2] <= 0.) return 0x0;
+
+ fdEdx[0] = TMath::Log(dedx[0]+dedx[1]);
+ fdEdx[1] = TMath::Log(dedx[2]);
+ return fdEdx;
}
-//__________________________________________________________________
-// Double_t AliTRDpidRefMakerLQ::Estimate3D2(TH3 * const h, Float_t &x, Float_t &y, Float_t &z)
-// {
-// // Author Alex Bercuci (A.Bercuci@gsi.de)
-// return 0.;
-// }
-
-
-///////////// Private functions ///////////////////////////////////
-
-//__________________________________________________________________
-Int_t AliTRDpidRefMakerLQ::CheckProdDirTree(const Char_t *dir)
+//________________________________________________________________________
+void AliTRDpidRefMakerLQ::Fill()
{
- // Scan directory tree for momenta. Returns the smallest number of
- // batches found in all directories or 0 if one momentum is missing.
-
- const char *pwd = gSystem->pwd();
-
- if(!gSystem->ChangeDirectory(dir)){
- AliError(Form("Couldn't access production root directory %s.", dir));
- return 0;
- }
-
- Int_t iDir;
- Int_t nDir = Int_t(1.E6);
- for(int imom=0; imom<AliTRDCalPID::kNMom; imom++){
- if(!gSystem->ChangeDirectory(Form("%3.1fGeV", AliTRDCalPID::GetMomentum(imom)))){
- AliError(Form("Couldn't find data for momentum %3.1f GeV/c.", AliTRDCalPID::GetMomentum(imom)));
- return 0;
- }
-
- iDir = 0;
- while(gSystem->ChangeDirectory(Form("%03d", iDir))){
- iDir++;
- gSystem->ChangeDirectory("..");
- }
- if(iDir < nDir) nDir = iDir;
- gSystem->ChangeDirectory(dir);
- }
-
- gSystem->ChangeDirectory(pwd);
-
- return nDir;
+ // get particle type
+ fSbin = TMath::LocMax(AliPID::kSPECIES, fPID);
+ // get momentum bin
+ fPbin = AliTRDpidUtil::GetMomentumBin(fP);
+ // fill data
+ fData->Fill();
}
-
-//__________________________________________________________________
-void AliTRDpidRefMakerLQ::Prepare2D()
+//________________________________________________________________________
+Bool_t AliTRDpidRefMakerLQ::PostProcess()
{
- //
- // Prepares the 2-dimensional reference histograms
- //
-
- // histogram settings
- Float_t xmin, xmax, ymin, ymax;
- Int_t nbinsx, nbinsy, nBinsSector, nSectors;
- const Int_t color[] = {4, 3, 2, 7, 6};
- for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
- // check PCA data
- if(!fPrinc[ispec]){
- AliError(Form("No data defined for %s.", AliTRDCalPID::GetPartName(ispec)));
- return;
- }
- // build reference histograms
- nBinsSector = 10;
- nSectors = 20;
- nbinsx = nbinsy = nBinsSector * nSectors;
- xmin = ymin = 0.;
- xmax = 8000.; ymax = 6000.;
- if(!fH2dEdx[ispec]){
- fH2dEdx[ispec] = new TH2D(Form("h2%s", AliTRDCalPID::GetPartSymb(ispec)), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
- fH2dEdx[ispec]->SetLineColor(color[ispec]);
- }
+ TFile::Open(Form("TRD.Calib%s.root", GetName()));
+ fData = dynamic_cast<TTree*>(gFile->Get(GetName()));
+ if (!fData) {
+ AliError("Tree not available");
+ return kFALSE;
}
+ AliDebug(2, Form("Data[%d]", fData->GetEntries()));
- // build transformed rotated histograms
- nbinsx = nbinsy = 100;
- xmin = ymin = -6.;
- xmax = 10.; ymax = 6.;
- TH2I *hProj = 0x0;
- if(!(hProj = (TH2I*)gDirectory->Get("hProj"))) hProj = new TH2I("hProj", "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
- else hProj->Reset();
-
- printf("Doing interpolation and invertion ... "); fflush(stdout);
- Bool_t kDebugPlot = kTRUE;
- //Bool_t kDebugPrint = kFALSE;
-
- TCanvas *c = 0x0;
- TEllipse *ellipse = 0x0;
- TMarker *mark = 0x0;
- if(kDebugPlot){
- c=new TCanvas("c2", "Interpolation 2D", 10, 10, 500, 500);
- ellipse = new TEllipse();
- ellipse->SetFillStyle(0); ellipse->SetLineColor(2);
- mark = new TMarker();
- mark->SetMarkerColor(2); mark->SetMarkerSize(2); mark->SetMarkerStyle(2);
- }
-
- // define observable variables
- TAxis *ax, *ay;
- Double_t xy[2], lxy[2], rxy[2];
- Double_t estimate, position;
- const TVectorD *eValues;
-
- // define radial sectors
- const Int_t nPhi = 36;
- const Float_t dPhi = TMath::TwoPi()/nPhi;
- //const Float_t dPhiRange = .1;
- Int_t nPoints[nPhi], nFitPoints, binStart, binStop;
- TLinearFitter refsFitter[nPhi], refsLongFitter(6, "1++x++y++x*x++y*y++x*y");
- Float_t fgFitterRange[nPhi];
- Bool_t kFitterStatus[nPhi];
- for(int iphi=0; iphi<nPhi; iphi++){
- refsFitter[iphi].SetDim(3);
- refsFitter[iphi].SetFormula("1++x++y");//++x*x++y*y++x*y");
- fgFitterRange[iphi] = .8;
- kFitterStatus[iphi] = kFALSE;
- }
- std::vector<UShort_t> storeX[nPhi], storeY[nPhi];
-
- // define working variables
- Float_t x0, y0, rx, ry;
- //Float_t rc, rmin, rmax, dr, dCT;
- Double_t phi, r;
- Int_t iPhi;
- Double_t entries;
- for(int ispec=0; ispec<5; ispec++){
- hProj->Reset();
- for(int iphi=0; iphi<nPhi; iphi++){
- nPoints[iphi] = 0;
- refsFitter[iphi].ClearPoints();
- fgFitterRange[iphi] = .8;
- kFitterStatus[iphi] = kFALSE;
- storeX[iphi].clear();
- storeY[iphi].clear();
- }
-
- // calculate covariance ellipse
- fPrinc[ispec]->MakePrincipals();
- eValues = fPrinc[ispec]->GetEigenValues();
- x0 = 0.;
- y0 = 0.;
- rx = 3.5*sqrt((*eValues)[0]);
- ry = 3.5*sqrt((*eValues)[1]);
-
- // rotate to principal axis
- Int_t irow = 0;
- const Double_t *xx;
- printf("filling for spec %d ...\n", ispec);
- while((xx = fPrinc[ispec]->GetRow(irow++))){
- fPrinc[ispec]->X2P(xx, rxy);
- hProj->Fill(rxy[0], rxy[1]);
- }
-
-
-// // debug plot
-// if(kDebugPlot){
-// hProj->Draw();
-// ellipse->DrawEllipse(x0, y0, rx, ry, 0., 360., 0.);
-// mark->DrawMarker(x0, y0);
-// gPad->Modified(); gPad->Update();
-// }
-
- // define radial sectors
- ax=hProj->GetXaxis();
- ay=hProj->GetYaxis();
- for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
- rxy[0] = ax->GetBinCenter(ibin);
- for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
- rxy[1] = ay->GetBinCenter(jbin);
-
- if((entries = hProj->GetBinContent(ibin, jbin)) == 0) continue;
-
- position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
- if(position < 1.) continue;
-
- r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]);
- phi = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi]
- iPhi = nPhi/2 + Int_t(phi/dPhi) - ((phi/dPhi > 0.) ? 0 : 1);
-
- refsFitter[iPhi].AddPoint(rxy, log(entries), 1./sqrt(entries));
- nPoints[iPhi]++;
- }
- }
-
- // do interpolation
- ax=fH2dEdx[ispec]->GetXaxis();
- ay=fH2dEdx[ispec]->GetYaxis();
- for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
- xy[0] = ax->GetBinCenter(ibin);
- lxy[0] = log(xy[0]);
- for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
- xy[1] = ay->GetBinCenter(jbin);
- lxy[1] = log(xy[1]);
-
- // rotate to PCA
- fPrinc[ispec]->X2P(lxy, rxy);
-
- // calculate border of covariance ellipse
- position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
-
- // interpolation inside the covariance ellipse
- if(position < 1.){
- Float_t xTemp = rxy[0];
- Float_t yTemp = rxy[1];
- estimate = Estimate2D2((TH2 *) hProj, xTemp, yTemp);
- rxy[0] = xTemp;
- rxy[1] = yTemp;
- // estimate = Estimate2D2((TH2 *) hProj, (Float_t)rxy[0], (Float_t)rxy[1]);
- fH2dEdx[ispec]->SetBinContent(ibin, jbin, estimate/xy[0]/xy[1]);
- } else { // interpolation outside the covariance ellipse
- r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]);
- phi = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi]
- iPhi = nPhi/2 + Int_t(phi/dPhi) - ((phi/dPhi > 0.) ? 0 : 1);
-
- storeX[iPhi].push_back(ibin);
- storeY[iPhi].push_back(jbin);
- }
- }
- }
-
- // Fill outliers
- // Radial fit on transformed rotated
- TVectorD vec(3);
- Int_t xbin, ybin;
- for(int iphi=0; iphi<nPhi; iphi++){
- phi = iphi * dPhi - TMath::Pi();
- if(TMath::Abs(TMath::Abs(phi)-TMath::Pi()) < 100.*TMath::DegToRad()) continue;
-
-
- refsFitter[iphi].Eval();
- refsFitter[iphi].GetParameters(vec);
- for(UInt_t ipoint=0; ipoint<storeX[iphi].size(); ipoint++){
- xbin = storeX[iphi].at(ipoint);
- ybin = storeY[iphi].at(ipoint);
- xy[0] = ax->GetBinCenter(xbin); lxy[0] = log(xy[0]);
- xy[1] = ay->GetBinCenter(ybin); lxy[1] = log(xy[1]);
-
- // rotate to PCA
- fPrinc[ispec]->X2P(lxy, rxy);
- estimate = exp(vec[0] + rxy[0]*vec[1] + rxy[1]*vec[2]);//+ rxy[0]*rxy[0]*vec[3]+ rxy[1]*rxy[1]*vec[4]+ rxy[0]*rxy[1]*vec[5]);
- fH2dEdx[ispec]->SetBinContent(xbin, ybin, estimate/xy[0]/xy[1]);
- }
- }
-
- // Longitudinal fit on ref histo in y direction
- for(int isector=1; isector<nSectors; isector++){
- // define sectors
- binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
- binStop = (isector+1)*nBinsSector;
-
- nPoints[0] = 0;
- refsLongFitter.ClearPoints();
- storeX[0].clear();
- storeY[0].clear();
-
- for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
- xy[0] = ax->GetBinCenter(ibin);
- for(int jbin=binStart; jbin<=binStop; jbin++){
- xy[1] = ay->GetBinCenter(jbin);
- if((entries = fH2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
- refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
- nPoints[0]++;
- } else {
- storeX[0].push_back(ibin);
- storeY[0].push_back(jbin);
- }
- }
- if(nPoints[0] >= ((isector==0)?60:420)) break;
- }
-
-
- if(refsLongFitter.Eval() == 1){
- printf("Error in Y sector %d\n", isector);
+ Float_t *data[] = {0x0, 0x0};
+ TPrincipal principal(2, "ND");
+
+ TCanvas *cc = new TCanvas("cc", "", 500, 500);
+
+ for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
+ for(Int_t is=AliPID::kSPECIES; is--;) {
+ principal.Clear();
+ Int_t n = fData->Draw("dEdx[0]:dEdx[1]", Form("p==%d&&s==%d", ip, is), "goff");
+ AliDebug(2, Form("pBin[%d] sBin[%d] n[%d]", ip, is, n));
+ if(n<10){
+ AliWarning(Form("Not enough data for %s[%d].", AliPID::ParticleShortName(is), ip));
continue;
}
- refsLongFitter.GetParameters(vec);
- for(UInt_t ipoint=0; ipoint<storeX[0].size(); ipoint++){
- xbin = storeX[0].at(ipoint);
- ybin = storeY[0].at(ipoint);
- xy[0] = ax->GetBinCenter(xbin);
- xy[1] = ay->GetBinCenter(ybin);
-
- estimate = vec[0] + xy[0]*vec[1] + xy[1]*vec[2] + xy[0]*xy[0]*vec[3]+ xy[1]*xy[1]*vec[4]+ xy[0]*xy[1]*vec[5];
- fH2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
+ // allocate storage
+ data[0] = new Float_t[n];data[1] = new Float_t[n];
+
+ // fill and make principal
+ Double_t *v1 = fData->GetV1(), *v2 = fData->GetV2();
+ while(n--){
+ Double_t dedx[] = {v1[n], v2[n]};
+ principal.AddRow(dedx);
+ }
+ principal.MakePrincipals();
+
+// // calculate covariance ellipse
+// eValues = principal.GetEigenValues();
+// x0 = 0.;
+// y0 = 0.;
+// rx = 3.5*sqrt((*eValues)[0]);
+// ry = 3.5*sqrt((*eValues)[1]);
+
+ // rotate to principal axis
+ const Double_t *xx; Double_t rxy[2];
+ while((xx = principal.GetRow(++n))){
+ principal.X2P(xx, rxy);
+ data[0][n]=rxy[0]; data[1][n]=rxy[1];
+ //hProj->Fill(rxy[0], rxy[1]);
}
- }
- // Longitudinal fit on ref histo in x direction
- for(int isector=1; isector<nSectors; isector++){
- binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
- binStop = (isector+1)*nBinsSector;
-
- nPoints[0] = 0;
- nFitPoints = 0;
- refsLongFitter.ClearPoints();
- storeX[0].clear();
- storeY[0].clear();
-
- for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
- xy[1] = ay->GetBinCenter(jbin);
- for(int ibin=binStart; ibin<=binStop; ibin++){
- xy[0] = ax->GetBinCenter(ibin);
- if((entries = fH2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
- refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
- nPoints[0]++;
- } else {
- storeX[0].push_back(ibin);
- storeY[0].push_back(jbin);
- nFitPoints++;
- }
- }
- if(nPoints[0] >= ((isector==0)?60:420)) break;
+ // estimate acceptable statistical error per bucket
+
+ // estimate number of buckets
+
+ // build PDF
+ TKDPDF pdf(n, 2, 10, data);
+ printf("PDF nodes[%d]\n", pdf.GetNTNodes());
+ pdf.SetStore();
+ pdf.SetAlpha(5.);
+ Float_t *c, v, ve; Double_t r, e;
+ for(Int_t in=pdf.GetNTNodes(); in--;){
+ pdf.GetCOGPoint(in, c, v, ve);
+ rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
+ printf("%2d x[%f] y[%f]\n", in, rxy[0], rxy[1]);
+ pdf.Eval(rxy, r, e, kTRUE);
}
- if(nFitPoints == 0) continue;
-
+ pdf.DrawBins(0,1,-6,6,-6,6);cc->Modified(); cc->Update();
+ AliDebug(2, Form("n[%d]", n));
+
+ //pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
- if(refsLongFitter.Eval() == 1){
- printf("Error in X sector %d\n", isector);
- continue;
- }
- refsLongFitter.GetParameters(vec);
- for(UInt_t ipoint=0; ipoint<storeX[0].size(); ipoint++){
- xbin = storeX[0].at(ipoint);
- ybin = storeY[0].at(ipoint);
- xy[0] = ax->GetBinCenter(xbin);
- xy[1] = ay->GetBinCenter(ybin);
-
- estimate = vec[0] + xy[0]*vec[1] + xy[1]*vec[2] + xy[0]*xy[0]*vec[3]+ xy[1]*xy[1]*vec[4]+ xy[0]*xy[1]*vec[5];
- fH2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
- }
- }
-
- fH2dEdx[ispec]->Scale(1./fH2dEdx[ispec]->Integral());
-
- // debug plot
- if(kDebugPlot){
- fH2dEdx[ispec]->Draw("cont1"); gPad->SetLogz();
- gPad->Modified(); gPad->Update();
+ delete [] data[0]; delete [] data[1];
}
}
- AliInfo("Finish interpolation.");
+
+ return kTRUE; // testing protection
}
+
//__________________________________________________________________
void AliTRDpidRefMakerLQ::Reset()
{
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
if(fH2dEdx[ispec]) fH2dEdx[ispec]->Reset();
- fPrinc[ispec]->Clear();
}
- if(fH1TB[0]) fH1TB[0]->Reset();
- if(fH1TB[1]) fH1TB[1]->Reset();
}
//__________________________________________________________________
if(!kFOUND) fSave = TFile::Open(fn, "RECREATE");
fSave->cd();
- Float_t fmom = AliTRDCalPID::GetMomentum(mom);
-
// save dE/dx references
TH2 *h2 = 0x0;
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
h2->Write();
}
- // save maximum time bin references
- TH1 *h1 = 0x0;
- h1 = (TH1F*)fH1TB[0]->Clone(Form("h1MaxTBEL%02d", mom));
- h1->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fmom));
- h1->GetXaxis()->SetTitle("time [100 ns]");
- h1->GetYaxis()->SetTitle("Probability");
- h1->Write();
-
- h1 = (TH1F*)fH1TB[1]->Clone(Form("h1MaxTBPI%02d", mom));
- h1->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fmom));
- h1->GetXaxis()->SetTitle("time [100 ns]");
- h1->GetYaxis()->SetTitle("Probability");
- h1->Write();
-
-
pwd->cd();
}
// //
///////////////////////////////////////////////////////////////////////////////
+#ifndef ALIPID_H
+#include "AliPID.h"
+#endif
+
#ifndef ALITRDPIDREFMAKER_H
#include "AliTRDpidRefMaker.h"
#endif
-class TH1;
-class TH2;
-class TH3;
-class TPrincipal;
-class TLinearFitter;
+#ifndef ALITRDPIDUTIL_H
+#include "AliTRDpidUtil.h"
+#endif
+class TH2;
class AliTRDpidRefMakerLQ : public AliTRDpidRefMaker {
-
public:
+ enum ETRDpidRefMakerLQsteer{
+ kMaxStat = 100000 // maximum statistics/species/momentum
+ };
AliTRDpidRefMakerLQ();
~AliTRDpidRefMakerLQ();
- static Double_t Estimate2D2(TH2 * const h, Float_t &x, Float_t &y);
- static Double_t Estimate2D1(TH2 * const h, Float_t &x, Float_t &y, const Float_t &dCT
- , const Float_t &rmin, const Float_t &rmax);
- // Double_t Estimate3D2(TH3 * const h, Float_t &x, Float_t &y, Float_t &z);
+ void CreateOutputObjects();
+ Bool_t PostProcess();
protected:
- void MakeRefs(Int_t pbin);
+ Float_t* GetdEdx(AliTRDseedV1*);
+ Int_t GetNslices() { return 2;}
+ void Fill();
private:
AliTRDpidRefMakerLQ(const AliTRDpidRefMakerLQ &ref);
AliTRDpidRefMakerLQ& operator=(const AliTRDpidRefMakerLQ &ref);
- Int_t CheckProdDirTree(const Char_t *dir=".");
- void Prepare2D();
- void Reset();
- void SaveReferences(const Int_t mom, const char *fn);
+ void Reset();
+ void SaveReferences(const Int_t mom, const char *fn);
private:
- TPrincipal *fPrinc[5]; // Used for principal component analysis
- static TLinearFitter *fgFitter2D2; // Working object for linear fitter
- static TLinearFitter *fgFitter2D1; // Working object for linear fitter
- TH2 *fH2dEdx[5]; // dE/dx data holders
- TH1 *fH1TB[2]; // Max time bin data holders
+ UChar_t fPbin; //! momentum bin
+ UChar_t fSbin; //! species bin
+ TH2 *fH2dEdx[5]; //! dE/dx data holders
- ClassDef(AliTRDpidRefMakerLQ, 3) // Reference histograms builder
+ ClassDef(AliTRDpidRefMakerLQ, 4) // Reference builder for Multidim-LQ TRD-PID
};
#include "TSystem.h"
#include "TDatime.h"
-#include "TEventList.h"
+#include "TPDGCode.h"
+#include "TH1F.h"
+#include "TFile.h"
#include "TGraphErrors.h"
+#include "TTree.h"
#include "TTreeStream.h"
-#include "TH1F.h"
+#include "TEventList.h"
#include "TMultiLayerPerceptron.h"
-#include "AliLog.h"
#include "AliPID.h"
#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliTrackReference.h"
+
+#include "AliAnalysisTask.h"
+#include "AliTRDtrackV1.h"
+#include "AliTRDReconstructor.h"
+#include "../Cal/AliTRDCalPID.h"
#include "../Cal/AliTRDCalPIDNN.h"
#include "AliTRDpidUtil.h"
#include "AliTRDpidRefMakerNN.h"
#include "info/AliTRDtrackInfo.h"
+#include "info/AliTRDv0Info.h"
// builds the reference tree for the training of neural networks
//________________________________________________________________________
AliTRDpidRefMakerNN::AliTRDpidRefMakerNN()
- :AliTRDpidRefMaker(UChar_t(AliTRDpidUtil::kNNslices), "PidRefMakerNN", "PID(NN) Reference Maker")
+ :AliTRDrecoTask("PidRefMakerNN", "PID(NN) Reference Maker")
+ ,fReconstructor(0x0)
+ ,fV0s(0x0)
+ ,fNN(0x0)
+ ,fLayer(0xff)
,fTrainMomBin(kAll)
,fEpochs(1000)
,fMinTrain(100)
,fDate(0)
+ ,fMom(0.)
,fDoTraining(0)
,fContinueTraining(0)
,fTrainPath(0x0)
//
// Default constructor
//
- SetAbundance(.67, .33);
- SetScaledEdx(Float_t(AliTRDCalPIDNN::kMLPscale));
+
+ fReconstructor = new AliTRDReconstructor();
+ fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
+ memset(fv0pid, 0, AliPID::kSPECIES*sizeof(Float_t));
+ memset(fdEdx, 0, 10*sizeof(Float_t));
+
+ const Int_t nnSize = AliTRDCalPID::kNMom * AliTRDgeometry::kNlayer;
+ memset(fTrain, 0, nnSize*sizeof(TEventList*));
+ memset(fTest, 0, nnSize*sizeof(TEventList*));
+ memset(fNet, 0, AliTRDgeometry::kNlayer*sizeof(TMultiLayerPerceptron*));
+
TDatime datime;
fDate = datime.GetDate();
+
+ DefineInput(1, TObjArray::Class());
+ DefineOutput(1, TTree::Class());
}
//________________________________________________________________________
AliTRDpidRefMakerNN::~AliTRDpidRefMakerNN()
{
+ if(fReconstructor) delete fReconstructor;
+ //if(fNN) delete fNN;
+ //if(fLQ) delete fLQ;
}
+//________________________________________________________________________
+void AliTRDpidRefMakerNN::ConnectInputData(Option_t *opt)
+{
+ AliTRDrecoTask::ConnectInputData(opt);
+ fV0s = dynamic_cast<TObjArray*>(GetInputData(1));
+}
+
//________________________________________________________________________
void AliTRDpidRefMakerNN::CreateOutputObjects()
{
// Create histograms
// Called once
- AliTRDpidRefMaker::CreateOutputObjects();
+
+ OpenFile(0, "RECREATE");
+ fContainer = new TObjArray();
+ fContainer->SetName(Form("Moni%s", GetName()));
+ fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5), kHistoPDG);
+
TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
+ gEffisTrain->SetNameTitle("EffTrain", "Train Efficiency Monitor");
gEffisTrain -> SetLineColor(4);
gEffisTrain -> SetMarkerColor(4);
gEffisTrain -> SetMarkerStyle(29);
gEffisTrain -> SetMarkerSize(1);
TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
+ gEffisTest->SetNameTitle("EffTest", "Test Efficiency Monitor");
gEffisTest -> SetLineColor(2);
gEffisTest -> SetMarkerColor(2);
gEffisTest -> SetMarkerStyle(29);
fContainer -> AddAt(gEffisTrain,kGraphTrain);
fContainer -> AddAt(gEffisTest,kGraphTest);
+
+ // open reference TTree for NN
+ fNN = new TTree("NN", "Reference data for NN");
+ fNN->Branch("fLayer", &fLayer, "fLayer/I");
+ fNN->Branch("fMom", &fMom, "fMom/F");
+ fNN->Branch("fv0pid", fv0pid, Form("fv0pid[%d]/F", AliPID::kSPECIES));
+ fNN->Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDpidUtil::kNNslices));
+}
+
+
+//________________________________________________________________________
+void AliTRDpidRefMakerNN::Exec(Option_t *)
+{
+ // Main loop
+ // Called for each event
+
+ Int_t labelsacc[10000];
+ memset(labelsacc, 0, sizeof(Int_t) * 10000);
+
+ Float_t mom;
+ ULong_t status;
+ Int_t nTRD = 0;
+
+ AliTRDtrackInfo *track = 0x0;
+ AliTRDv0Info *v0 = 0x0;
+ AliTRDtrackV1 *TRDtrack = 0x0;
+ AliTrackReference *ref = 0x0;
+ AliExternalTrackParam *esd = 0x0;
+ AliTRDseedV1 *TRDtracklet = 0x0;
+
+ for(Int_t iv0=0; iv0<fV0s->GetEntriesFast(); iv0++){
+ v0 = dynamic_cast<AliTRDv0Info*>(fV0s->At(iv0));
+ v0->Print();
+ }
+
+ for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
+
+ // reset the pid information
+ for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++)
+ fv0pid[iPart] = 0.;
+
+ track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
+ if(!track->HasESDtrack()) continue;
+ status = track->GetStatus();
+ if(!(status&AliESDtrack::kTPCout)) continue;
+
+ if(!(TRDtrack = track->GetTrack())) continue;
+ //&&(track->GetNumberOfClustersRefit()
+
+ // use only tracks that hit 6 chambers
+ if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
+
+ ref = track->GetTrackRef(0);
+ esd = track->GetESDinfo()->GetOuterParam();
+ mom = ref ? ref->P(): esd->P();
+ fMom = mom;
+
+
+ labelsacc[nTRD] = track->GetLabel();
+ nTRD++;
+
+ // if no monte carlo data available -> use V0 information
+ if(!HasMCdata()){
+ GetV0info(TRDtrack,fv0pid);
+ }
+ // else use the MC info
+ else{
+ switch(track -> GetPDG()){
+ case kElectron:
+ case kPositron:
+ fv0pid[AliPID::kElectron] = 1.;
+ break;
+ case kMuonPlus:
+ case kMuonMinus:
+ fv0pid[AliPID::kMuon] = 1.;
+ break;
+ case kPiPlus:
+ case kPiMinus:
+ fv0pid[AliPID::kPion] = 1.;
+ break;
+ case kKPlus:
+ case kKMinus:
+ fv0pid[AliPID::kKaon] = 1.;
+ break;
+ case kProton:
+ case kProtonBar:
+ fv0pid[AliPID::kProton] = 1.;
+ break;
+ }
+ }
+
+ // set reconstructor
+ Float_t *dedx;
+ TRDtrack -> SetReconstructor(fReconstructor);
+
+ // fill the dE/dx information for NN
+ fReconstructor -> SetOption("nn");
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
+ if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
+ TRDtracklet->CookdEdx(AliTRDpidUtil::kNNslices);
+ dedx = const_cast<Float_t *>(TRDtracklet->GetdEdx());
+ for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kNNslices; iSlice++)
+ dedx[iSlice] = dedx[iSlice]/AliTRDCalPIDNN::kMLPscale;
+ memcpy(fdEdx, dedx, AliTRDpidUtil::kNNslices*sizeof(Float_t));
+ if(fDebugLevel>=2) Printf("LayerNN : %d", ily);
+ fLayer = ily;
+ fNN->Fill();
+ }
+
+
+
+ for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
+ if(fDebugLevel>=4) Printf("PDG is %d %f", iPart, fv0pid[iPart]);
+ }
+ }
+
+ PostData(0, fContainer);
+ PostData(1, fNN);
}
// Called once at the end of the query
// build the training andthe test list for the neural networks
- MakeTrainingList();
+ MakeTrainingLists();
if(!fDoTraining) return kTRUE;
// train the neural networks and build the refrence histos for 2-dim LQ
gSystem->Exec(Form("mkdir ./Networks_%d/",fDate));
- AliDebug(3, Form("TrainMomBin [%d] [%d]", fTrainMomBin, kAll));
+ if(fDebugLevel>=2) Printf("TrainMomBin [%d] [%d]", fTrainMomBin, kAll);
// train single network for a single momentum (recommended)
if(!(fTrainMomBin == kAll)){
if(fTrain[fTrainMomBin][0] -> GetN() < fMinTrain){
- AliWarning("Not enough events for training available! Please check Data sample!");
+ if(fDebugLevel>=2) Printf("Warning in AliTRDpidRefMakerNN::PostProcess : Not enough events for training available! Please check Data sample!");
return kFALSE;
}
- MakeRefs(fTrainMomBin);
+ TrainNetworks(fTrainMomBin);
MonitorTraining(fTrainMomBin);
}
// train all momenta
else{
for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
if(fTrain[iMomBin][0] -> GetN() < fMinTrain){
- AliWarning(Form("Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin));
- continue;
+ if(fDebugLevel>=2) Printf("Warning in AliTRDpidRefMakerNN::PostProcess : Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin);
+ continue;
}
- MakeRefs(fTrainMomBin);
+ TrainNetworks(iMomBin);
MonitorTraining(iMomBin);
}
}
}
+//________________________________________________________________________
+void AliTRDpidRefMakerNN::Terminate(Option_t *)
+{
+ // Draw result to the screen
+ // Called once at the end of the query
+
+ fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
+ if (!fContainer) {
+ Printf("ERROR: list not available");
+ return;
+ }
+}
+
+
+//________________________________________________________________________
+void AliTRDpidRefMakerNN::GetV0info(AliTRDtrackV1 *TRDtrack, Float_t *v0pid)
+{
+ // !!!! PREMILMINARY FUNCTION !!!!
+ //
+ // this is the place for the V0 procedure
+ // as long as there is no one implemented,
+ // just the probabilities
+ // of the TRDtrack are used!
+
+ TRDtrack -> SetReconstructor(fReconstructor);
+ fReconstructor -> SetOption("nn");
+ TRDtrack -> CookPID();
+ for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
+ v0pid[iPart] = TRDtrack -> GetPID(iPart);
+ if(fDebugLevel>=4) Printf("PDG is (in V0info) %d %f", iPart, v0pid[iPart]);
+ }
+}
+
//________________________________________________________________________
-void AliTRDpidRefMakerNN::MakeRefs(Int_t mombin)
+void AliTRDpidRefMakerNN::MakeTrainingLists()
+{
+ //
+ // build the training lists for the neural networks
+ //
+
+ if (!fNN) {
+ LoadFile("TRD.CalibPidRefMakerNN.root");
+ }
+
+ if (!fNN) {
+ Printf("ERROR tree for training list not available");
+ return;
+ }
+
+ if(fDebugLevel>=2) Printf("\n Making training lists! \n");
+
+ Int_t nPart[AliPID::kSPECIES][AliTRDCalPID::kNMom];
+ memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
+
+ // set needed branches
+ fNN -> SetBranchAddress("fv0pid", &fv0pid);
+ fNN -> SetBranchAddress("fMom", &fMom);
+ fNN -> SetBranchAddress("fLayer", &fLayer);
+
+ AliTRDpidUtil *util = new AliTRDpidUtil();
+
+ // start first loop to check total number of each particle type
+ for(Int_t iEv=0; iEv < fNN -> GetEntries(); iEv++){
+ fNN -> GetEntry(iEv);
+
+ // use only events with goes through 6 layers TRD
+ if(!fLayer == 0)
+ continue;
+
+ // set the 11 momentum bins
+ Int_t iMomBin = -1;
+ iMomBin = util -> GetMomentumBin(fMom);
+
+ // check PID information and count particle types per momentum interval
+ if(fv0pid[AliPID::kElectron] == 1)
+ nPart[AliPID::kElectron][iMomBin]++;
+ else if(fv0pid[AliPID::kMuon] == 1)
+ nPart[AliPID::kMuon][iMomBin]++;
+ else if(fv0pid[AliPID::kPion] == 1)
+ nPart[AliPID::kPion][iMomBin]++;
+ else if(fv0pid[AliPID::kKaon] == 1)
+ nPart[AliPID::kKaon][iMomBin]++;
+ else if(fv0pid[AliPID::kProton] == 1)
+ nPart[AliPID::kProton][iMomBin]++;
+ }
+
+ if(fDebugLevel>=2){
+ Printf("Particle multiplicities:");
+ for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++)
+ Printf("Momentum[%d] Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", iMomBin, nPart[AliPID::kElectron][iMomBin], nPart[AliPID::kMuon][iMomBin], nPart[AliPID::kPion][iMomBin], nPart[AliPID::kKaon][iMomBin], nPart[AliPID::kProton][iMomBin]);
+ Printf("\n");
+ }
+
+ // implement counter of training and test sample size
+ Int_t iTrain[AliTRDCalPID::kNMom], iTest[AliTRDCalPID::kNMom];
+ memset(iTrain, 0, AliTRDCalPID::kNMom*sizeof(Int_t));
+ memset(iTest, 0, AliTRDCalPID::kNMom*sizeof(Int_t));
+
+ // set training sample size per momentum interval to 2/3
+ // of smallest particle counter and test sample to 1/3
+ for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
+ iTrain[iMomBin] = nPart[0][iMomBin];
+ for(Int_t iPart = 1; iPart < AliPID::kSPECIES; iPart++){
+ if(iTrain[iMomBin] > nPart[iPart][iMomBin])
+ iTrain[iMomBin] = nPart[iPart][iMomBin];
+ }
+ iTrain[iMomBin] = Int_t(iTrain[iMomBin] * .66);
+ iTest[iMomBin] = Int_t( iTrain[iMomBin] * .5);
+ if(fDebugLevel>=2) Printf("Momentum[%d] Train[%d] Test[%d]", iMomBin, iTrain[iMomBin], iTest[iMomBin]);
+ }
+ if(fDebugLevel>=2) Printf("\n");
+
+
+ // reset couters
+ memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
+
+ // start second loop to set the event lists
+ for(Int_t iEv = 0; iEv < fNN -> GetEntries(); iEv++){
+ fNN -> GetEntry(iEv);
+
+ // use only events with goes through 6 layers TRD
+ if(!fLayer == 0)
+ continue;
+
+ // set the 11 momentum bins
+ Int_t iMomBin = -1;
+ iMomBin = util -> GetMomentumBin(fMom);
+
+ // set electrons
+ if(fv0pid[AliPID::kElectron] == 1){
+ if(nPart[AliPID::kElectron][iMomBin] < iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTrain[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kElectron][iMomBin]++;
+ }
+ else if(nPart[AliPID::kElectron][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTest[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kElectron][iMomBin]++;
+ }
+ else
+ continue;
+ }
+ // set muons
+ else if(fv0pid[AliPID::kMuon] == 1){
+ if(nPart[AliPID::kMuon][iMomBin] < iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTrain[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kMuon][iMomBin]++;
+ }
+ else if(nPart[AliPID::kMuon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTest[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kMuon][iMomBin]++;
+ }
+ else
+ continue;
+ }
+ // set pions
+ else if(fv0pid[AliPID::kPion] == 1){
+ if(nPart[AliPID::kPion][iMomBin] < iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTrain[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kPion][iMomBin]++;
+ }
+ else if(nPart[AliPID::kPion][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTest[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kPion][iMomBin]++;
+ }
+ else
+ continue;
+ }
+ // set kaons
+ else if(fv0pid[AliPID::kKaon] == 1){
+ if(nPart[AliPID::kKaon][iMomBin] < iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTrain[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kKaon][iMomBin]++;
+ }
+ else if(nPart[AliPID::kKaon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTest[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kKaon][iMomBin]++;
+ }
+ else
+ continue;
+ }
+ // set protons
+ else if(fv0pid[AliPID::kProton] == 1){
+ if(nPart[AliPID::kProton][iMomBin] < iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTrain[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kProton][iMomBin]++;
+ }
+ else if(nPart[AliPID::kProton][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
+ fTest[iMomBin][ily] -> Enter(iEv + ily);
+ nPart[AliPID::kProton][iMomBin]++;
+ }
+ else
+ continue;
+ }
+ }
+
+ if(fDebugLevel>=2){
+ Printf("Particle multiplicities in both lists:");
+ for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++)
+ Printf("Momentum[%d] Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", iMomBin, nPart[AliPID::kElectron][iMomBin], nPart[AliPID::kMuon][iMomBin], nPart[AliPID::kPion][iMomBin], nPart[AliPID::kKaon][iMomBin], nPart[AliPID::kProton][iMomBin]);
+ Printf("\n");
+ }
+
+ util -> Delete();
+}
+
+
+//________________________________________________________________________
+void AliTRDpidRefMakerNN::TrainNetworks(Int_t mombin)
{
//
// train the neural networks
//
- if (!fData) LoadFile(Form("TRD.%s.root", GetName()));
+ if (!fNN) {
+ LoadFile("TRD.CalibPidRefMakerNN.root");
+ }
- if (!fData) {
- AliError("Tree for training list not available");
+ if (!fNN) {
+ Printf("ERROR tree for training list not available");
return;
}
TDatime datime;
fDate = datime.GetDate();
- AliDebug(2, Form("Training momentum bin %d", mombin));
+ if(fDebugLevel>=2) Printf("Training momentum bin %d", mombin);
// set variable to monitor the training and to save the development of the networks
Int_t nEpochs = fEpochs/kMoniTrain;
- AliDebug(2, Form("Training %d times %d epochs", kMoniTrain, nEpochs));
+ if(fDebugLevel>=2) Printf("Training %d times %d epochs", kMoniTrain, nEpochs);
// make directories to save the networks
gSystem->Exec(Form("rm -r ./Networks_%d/MomBin_%d",fDate, mombin));
// loop over chambers
for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
// set the event lists
- fData -> SetEventList(fTrain[mombin][iChamb]);
- fData -> SetEventList(fTest[mombin][iChamb]);
+ fNN -> SetEventList(fTrain[mombin][iChamb]);
+ fNN -> SetEventList(fTest[mombin][iChamb]);
- AliDebug(2, Form("Trainingloop[%d] Chamber[%d]", iLoop, iChamb));
+ if(fDebugLevel>=2) Printf("Trainingloop[%d] Chamber[%d]", iLoop, iChamb);
// check if network is already implemented
if(bFirstLoop[iChamb] == kTRUE){
- fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fData,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
+ fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fv0pid[0],fv0pid[1],fv0pid[2],fv0pid[3],fv0pid[4]!",fNN,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
fNet[iChamb] -> SetLearningMethod(TMultiLayerPerceptron::kStochastic); // set learning method
fNet[iChamb] -> TMultiLayerPerceptron::SetEta(0.001); // set learning speed
if(!fContinueTraining){
}
+
//________________________________________________________________________
void AliTRDpidRefMakerNN::MonitorTraining(Int_t mombin)
{
// train the neural networks
//
- if(!fContainer) LoadContainer(Form("TRD.Task%s.root", GetName()));
-
if(!fContainer){
- AliError("ERROR container not available");
+ LoadContainer("TRD.CalibPidRefMakerNN.root");
+ }
+ if(!fContainer){
+ Printf("ERROR container not available");
return;
}
- if (!fData) LoadFile(Form("TRD.Task%s.root", GetName()));
- if (!fData) {
- AliError("Tree for training list not available");
+ if (!fNN) {
+ LoadFile("TRD.CalibPidRefMakerNN.root");
+ }
+ if (!fNN) {
+ Printf("ERROR tree for training list not available");
return;
}
// init networks and set event list
for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
- fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fData,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
- fData -> SetEventList(fTrain[mombin][iChamb]);
- fData -> SetEventList(fTest[mombin][iChamb]);
+ fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fv0pid[0],fv0pid[1],fv0pid[2],fv0pid[3],fv0pid[4]!",fNN,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
+ fNN -> SetEventList(fTrain[mombin][iChamb]);
+ fNN -> SetEventList(fTest[mombin][iChamb]);
}
// implement variables for likelihoods
}
TotProb = 0.;
- fData -> GetEntry(fTrain[mombin][0] -> GetEntry(iEvent));
+ fNN -> GetEntry(fTrain[mombin][0] -> GetEntry(iEvent));
// use event only if it is electron or pion
- if(!((fPID[AliPID::kElectron] == 1.0) || (fPID[AliPID::kPion] == 1.0))) continue;
+ if(!((fv0pid[AliPID::kElectron] == 1.0) || (fv0pid[AliPID::kPion] == 1.0))) continue;
// get the probabilities for each particle type in each chamber
for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
}
// fill likelihood distributions
- if(fPID[AliPID::kElectron] == 1)
+ if(fv0pid[AliPID::kElectron] == 1)
hElecs -> Fill(LikeAll[AliPID::kElectron]);
- if(fPID[AliPID::kPion] == 1)
+ if(fv0pid[AliPID::kPion] == 1)
hPions -> Fill(LikeAll[AliPID::kElectron]);
} // end event loop
}
TotProb = 0.;
- fData -> GetEntry(fTest[mombin][0] -> GetEntry(iEvent));
+ fNN -> GetEntry(fTest[mombin][0] -> GetEntry(iEvent));
// use event only if it is electron or pion
- if(!((fPID[AliPID::kElectron] == 1.0) || (fPID[AliPID::kPion] == 1.0))) continue;
+ if(!((fv0pid[AliPID::kElectron] == 1.0) || (fv0pid[AliPID::kPion] == 1.0))) continue;
// get the probabilities for each particle type in each chamber
for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
}
// fill likelihood distributions
- if(fPID[AliPID::kElectron] == 1)
+ if(fv0pid[AliPID::kElectron] == 1)
hElecs -> Fill(LikeAll[AliPID::kElectron]);
- if(fPID[AliPID::kPion] == 1)
+ if(fv0pid[AliPID::kPion] == 1)
hPions -> Fill(LikeAll[AliPID::kElectron]);
} // end event loop
}
+//________________________________________________________________________
+void AliTRDpidRefMakerNN::LoadFile(const Char_t *InFileNN)
+{
+ //
+ // Loads the files and sets the event list
+ // for neural network training and
+ // building of the 2-dim reference histograms.
+ // Useable for training outside of the makeResults.C macro
+ //
+
+ TFile *fInFileNN;
+ fInFileNN = new TFile(InFileNN, "READ");
+ fNN = (TTree*)fInFileNN -> Get("NN");
+
+ for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
+ for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
+ fTrain[iMom][ily] = new TEventList(Form("fTrainMom%d_%d", iMom, ily), Form("Training list for momentum intervall %d and plane %d", iMom, ily));
+ fTest[iMom][ily] = new TEventList(Form("fTestMom%d_%d", iMom, ily), Form("Test list for momentum intervall %d and plane %d", iMom, ily));
+ }
+ }
+}
+
+
+//________________________________________________________________________
+void AliTRDpidRefMakerNN::LoadContainer(const Char_t *InFileCont)
+{
+
+ //
+ // Loads the container if no container is there.
+ // Useable for training outside of the makeResults.C macro
+ //
+
+ TFile *fInFileCont;
+ fInFileCont = new TFile(InFileCont, "READ");
+ fContainer = (TObjArray*)fInFileCont -> Get("PidRefMaker");
+
+}
+
+
+// //________________________________________________________________________
+// void AliTRDpidRefMakerNN::CreateGraphs()
+// {
+// // Create histograms
+// // Called once
+
+// OpenFile(0, "RECREATE");
+// fContainer = new TObjArray();
+// fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
+
+// TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
+// gEffisTrain -> SetLineColor(4);
+// gEffisTrain -> SetMarkerColor(4);
+// gEffisTrain -> SetMarkerStyle(29);
+// gEffisTrain -> SetMarkerSize(2);
+
+// TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
+// gEffisTest -> SetLineColor(2);
+// gEffisTest -> SetMarkerColor(2);
+// gEffisTest -> SetMarkerSize(2);
+
+// fContainer -> AddAt(gEffisTrain,kGraphTrain);
+// fContainer -> AddAt(gEffisTest,kGraphTest);
+// }
+
//
///////////////////////////////////////////////////////
-#ifndef ALITRDPIDREFMAKER_H
-#include "AliTRDpidRefMaker.h"
+#ifndef ALITRDRECOTASK_H
+#include "AliTRDrecoTask.h"
#endif
+#ifndef ALIPID_H
+#include "AliPID.h"
+#endif
+
+#ifndef ALITRDCALPID_H
+#include "../Cal/AliTRDCalPID.h"
+#endif
+
+#ifndef ALITRDGEOMETRY_H
+#include "../AliTRDgeometry.h"
+#endif
+
+class TTree;
+class TObjArray;
+class TEventList;
class TMultiLayerPerceptron;
-class AliTRDpidRefMakerNN : public AliTRDpidRefMaker
+class AliPID;
+class AliTRDtrackV1;
+class AliTRDReconstructor;
+class AliTRDpidRefMakerNN : public AliTRDrecoTask
{
+
public:
- enum ETRDpidRefMakerNNgraph {
- kGraphTrain = 0
- ,kGraphTest = 1
+ enum {
+ k006 = 0
+ ,k008 = 1
+ ,k010 = 2
+ ,k015 = 3
+ ,k020 = 4
+ ,k030 = 5
+ ,k040 = 6
+ ,k050 = 7
+ ,k060 = 8
+ ,k080 = 9
+ ,k100 = 10
+ ,kAll = 11
+ };
+
+ enum {
+ kHistoPDG = 0
+ ,kGraphTrain = 1
+ ,kGraphTest = 2
};
- enum ETRDpidRefMakerNNmoni {
+ enum {
kMoniTrain = 50
};
virtual ~AliTRDpidRefMakerNN();
+ void ConnectInputData(Option_t *opt);
void CreateOutputObjects();
-
+ void Exec(Option_t *option);
Int_t GetEpochs() {return fEpochs;};
Int_t GetMinTrain() {return fMinTrain;};
Int_t GetTrainMomBin() {return fTrainMomBin;};
void SetDoTraining(Bool_t train) {fDoTraining = train;};
void SetContinueTraining(Bool_t continTrain) {fContinueTraining = continTrain;};
void SetTrainPath(Int_t path) {fTrainPath = path;};
+ void LoadFile(const Char_t *InFileNN);
+ void Terminate(Option_t *);
+ void MakeTrainingLists(); // build the training and the test list
void MonitorTraining(Int_t mombin); // monitor training process
-
-protected:
- void MakeRefs(Int_t mombin); // train the neural networks for a given momentum bin
+ void LoadContainer(const Char_t *InFileCont);
+ //void CreateGraphs();
private:
AliTRDpidRefMakerNN(const AliTRDpidRefMakerNN&); // not implemented
AliTRDpidRefMakerNN& operator=(const AliTRDpidRefMakerNN&); // not implemented
+ void GetV0info(AliTRDtrackV1 *TRDtrack, Float_t *v0pdg); // get the v0 information
+ void TrainNetworks(Int_t mombin); // train the neural networks for a given momentum bin
+ AliTRDReconstructor *fReconstructor; //! reconstructor needed for recalculation the PID
+ TObjArray *fV0s; //! v0 array
+ TTree *fNN; // NN data
+ TEventList *fTrain[AliTRDCalPID::kNMom][AliTRDgeometry::kNlayer]; // Training list for each momentum
+ TEventList *fTest[AliTRDCalPID::kNMom][AliTRDgeometry::kNlayer]; // Test list for each momentum
TMultiLayerPerceptron *fNet[AliTRDgeometry::kNlayer]; // artificial neural network
+ Int_t fLayer; // TRD layer index
Int_t fTrainMomBin; // momentum bin for the training
Int_t fEpochs; // Number of epochs for the training of the NNs
Int_t fMinTrain; // minimum of events needed for training
Int_t fDate; // date stamp for training of the NNs
+ Float_t fMom; // momentum
+ Float_t *fdEdx[10]; // dEdx array
+ Float_t fv0pid[AliPID::kSPECIES]; // pid from v0s
Bool_t fDoTraining; // checks if training will be done
Bool_t fContinueTraining; // checks if training from an older run should be continued
Int_t fTrainPath; // sets the path for continuing the training
- ClassDef(AliTRDpidRefMakerNN, 2); // TRD PID reference maker for NN
+ ClassDef(AliTRDpidRefMakerNN, 2); // TRD reference maker for NN
};
#endif
if(!(TSTBIT(map, kPIDRefMaker))) return;
// TRD pid reference
- AliTRDpidRefMaker *ref = 0x0;
+ AliTRDpidRefMakerNN *ref = new AliTRDpidRefMakerNN();
+ //AliLog::SetClassDebugLevel("AliTRDpidRefMakerNN", 4);
// Neural network PID
- mgr->AddTask(ref = new AliTRDpidRefMakerNN());
+ mgr->AddTask(ref);
ref->SetDebugLevel(3);
AliLog::SetClassDebugLevel("AliTRDpidRefMakerNN", 3);
ref->SetMCdata(mgr->GetMCtruthEventHandler());
mgr->ConnectInput( ref, 0, ci[0]);
mgr->ConnectInput( ref, 1, ci[2]);
- mgr->ConnectOutput(ref, 0, mgr->CreateContainer(ref->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%s.root", ref->GetName())));
- mgr->ConnectOutput(ref, 1, mgr->CreateContainer(Form("%sNN", ref->GetName()), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.%s.root", ref->GetName())));
+ mgr->ConnectOutput(ref, 0, mgr->CreateContainer(Form("Moni%s", ref->GetName()), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", ref->GetName())));
+ mgr->ConnectOutput(ref, 1, mgr->CreateContainer(ref->GetName(), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", ref->GetName())));
// Multidimensional Likelihood PID
- mgr->AddTask(ref = new AliTRDpidRefMakerLQ());
- ref->SetDebugLevel(3);
- AliLog::SetClassDebugLevel("AliTRDpidRefMakerLQ", 3);
- ref->SetMCdata(mgr->GetMCtruthEventHandler());
- mgr->ConnectInput( ref, 0, ci[0]);
- mgr->ConnectInput( ref, 1, ci[2]);
- mgr->ConnectOutput(ref, 0, mgr->CreateContainer(ref->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%s.root", ref->GetName())));
- mgr->ConnectOutput(ref, 1, mgr->CreateContainer(Form("%sLQ", ref->GetName()), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.%s.root", ref->GetName())));
+ AliTRDpidRefMakerLQ *reflq = new AliTRDpidRefMakerLQ();
+ mgr->AddTask(reflq);
+ reflq->SetDebugLevel(3);
+ //AliLog::SetClassDebugLevel("AliTRDpidRefMakerLQ", 3);
+ reflq->SetMCdata(mgr->GetMCtruthEventHandler());
+ mgr->ConnectInput( reflq, 0, ci[0]);
+ mgr->ConnectInput( reflq, 1, ci[2]);
+ mgr->ConnectOutput(reflq, 0, mgr->CreateContainer(Form("Moni%s", reflq->GetName()), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", reflq->GetName())));
+ mgr->ConnectOutput(reflq, 1, mgr->CreateContainer(reflq->GetName(), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", reflq->GetName())));
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff