/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercialf purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ ///////////////////////////////////////////////////// // // Monitor V0 for TRD // // Authors: // Markus Heide ////////////////////////////////////////////////////// #include "TObjArray.h" #include "TH2.h" #include "TH2F.h" #include "TH1I.h" #include "TCanvas.h" #include "AliLog.h" #include "AliESDtrack.h" #include "AliPID.h" #include "AliTRDv0Monitor.h" #include "info/AliTRDv0Info.h" #include "info/AliTRDeventInfo.h" ClassImp(AliTRDv0Monitor) //________________________________________________________________________ AliTRDv0Monitor::AliTRDv0Monitor() :AliTRDrecoTask() ,fhQualityReductions(NULL) ,fV0s(NULL) ,fData(NULL) ,fInfo(NULL) ,fP(-1.) { // // Default constructor // SetNameTitle("TRDv0Monitor", "V0 Monitor for TRD PID"); } //________________________________________________________________________ AliTRDv0Monitor::AliTRDv0Monitor(const char *name) :AliTRDrecoTask(name, "V0 Monitor for TRD PID") ,fhQualityReductions(NULL) ,fV0s(NULL) ,fData(NULL) ,fInfo(NULL) ,fP(-1.) { // // Default constructor // DefineInput(3, TObjArray::Class()); // v0 list DefineInput(4, TObjArray::Class()); // pid info list } // //____________________________________________________________________ // void AliTRDv0Monitor::MakeSummary(){//makes a summary with potentially nice reference figures // //TCanvas *cOut = new TCanvas("v0MonitorSummary1", "Summary 1 for task V0Monitor", 1024, 768); // //cOut->cd(); // //GetRefFigure(4); // //cOut->SaveAs("V0MonitorSummary.gif"); // // //cOut = new TCanvas("v0MonitorSummary2","Summary 2 for task V0Monitor", 1024, 768); // //cOut->cd(); // //GetRefFigure(5); // //cOut->SaveAs("V0MonitorSummary2.gif"); // } //________________________________________________________________________ Bool_t AliTRDv0Monitor::GetRefFigure(Int_t /*ifig*/) { //creating reference figures AliInfo("Implementation on going ..."); return kTRUE; } //________________________________________________________________________ TObjArray* AliTRDv0Monitor::Histos() { // Create histograms // Called once if(fContainer) return fContainer; fContainer = new TObjArray(kNPlots); fContainer->SetOwner(); fContainer->SetName("V0Monitoring"); const char *samplename[AliPID::kSPECIES] = {"electrons","muons","pions","kaons","protons"}; const char *decayname[AliTRDv0Info::kNDecays] = {"gamma","K0s","Lambda","AntiLambda"}; const char *detectorname[kNDets] = {"ITS","TPC","TOF"}; fhQualityReductions = new TH1I(Form("fhQualityReductions"),Form("Number of tracks cut out by different quality cut steps"),11,-9,2); fContainer->Add(fhQualityReductions); for(Int_t ipart = 0;ipart < AliPID::kSPECIES; ipart++){ fhCutReductions[ipart] = new TH1I(Form("fhCutReductions_%s",samplename[ipart]),Form("Number of tracks cut out by different cut steps for %s",samplename[ipart]),19,-17,2); fContainer->Add(fhCutReductions[ipart]); for(Int_t idetector = 0; idetector < kNDets; idetector++){ fhDetPID[idetector][ipart] = new TH2F(Form("fhDetector_%s_%s",detectorname[idetector],samplename[ipart]),Form("%s Likelihood for %s vs. momentum",detectorname[idetector], samplename[ipart]),100,0.,10.,100, 0., 1.); fContainer->Add(fhDetPID[idetector][ipart]); } fhComPID[ipart] = new TH2F(Form("fhComPID_%s",samplename[ipart]),Form("Combined TPC/TOF PID: Likelihood for %s",samplename[ipart]),100,0.,10.,100,0.,1.); fContainer->Add(fhComPID[ipart]); for(Int_t cutstep = 0; cutstep < kNCutSteps; cutstep++){ fhTPCdEdx[ipart][cutstep] = new TH2F(Form("fhTPCdEdx_%s_[%d]",samplename[ipart],cutstep),Form("TPC dE/dx for %s [%d]",samplename[ipart],cutstep),100,0.,10.,300,0.,300.); fContainer->Add(fhTPCdEdx[ipart][cutstep]); } } for(Int_t iDecay = 0; iDecay < AliTRDv0Info::kNDecays; iDecay++){ for(Int_t cutstep =0; cutstep < kNCutSteps; cutstep++){ fhV0Chi2ndf[iDecay][cutstep] = new TH2F(Form("fhV0Chi2ndf_%s_[%d]",decayname[iDecay],cutstep),Form("Chi2/NDF vs. momentum"),100,0.,10.,500, 0., 500.); fContainer->Add(fhV0Chi2ndf[iDecay][cutstep]); fhPsiPair[iDecay][cutstep] = new TH2F(Form("fhV0PsiPair_%s_[%d]",decayname[iDecay],cutstep),Form("Psi_pair vs. momentum"),100,0.,10.,200, 0., 1.6); fContainer->Add(fhPsiPair[iDecay][cutstep]); fhPointAngle[iDecay][cutstep] = new TH2F(Form("fhPointAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Pointing Angle vs. momentum"),100,0.,10.,500, 0., 1.6); fContainer->Add(fhPointAngle[iDecay][cutstep]); fhDCA[iDecay][cutstep] = new TH2F(Form("fhDCA_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Daughter DCA vs. momentum"),100,0.,10.,500, 0., 1.); fContainer->Add(fhDCA[iDecay][cutstep]); fhOpenAngle[iDecay][cutstep] = new TH2F(Form("fhOpenAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Opening Angle vs. momentum"),100,0.,10.,500, 0., 1.6); fContainer->Add(fhOpenAngle[iDecay][cutstep]); fhRadius[iDecay][cutstep] = new TH2F(Form("fhRadius_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Generation Radius vs. momentum"),100,0.,10.,500, 0., 150.); fContainer->Add(fhRadius[iDecay][cutstep]); } fhInvMass[iDecay] = new TH2F(Form("fhInvMass_%s",decayname[iDecay]),Form("Invariant Mass vs. momentum"),100,0.,10.,500, 0., 2.); fContainer->Add(fhInvMass[iDecay]); } /*TH1F *hV0mcPID[AliPID::kSPECIES][AliPID::kSPECIES]; Int_t nPBins = 200; for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++){ for(Int_t iSample = 0; iSample < AliPID::kSPECIES; iSample++){ fhV0mcPID[iSample][iSpecies] = new TH1F(Form("fhV0mcPID_%s_is_%s",name[iSample],name[iSpecies]),Form("%s contained in %s sample",name[iSpecies],name[iSample]), nPBins, 0.2, 13.); } }*/ return fContainer; } //________________________________________________________________________ void AliTRDv0Monitor::UserExec(Option_t *) { // Main loop // Called for each event if(!(fTracks = dynamic_cast(GetInputData(1)))) return; if(!(fEvent = dynamic_cast(GetInputData(2)))) return; if(!(fV0s = dynamic_cast(GetInputData(3)))) return; if(!(fInfo = dynamic_cast(GetInputData(4)))) return; AliTRDtrackInfo *track = NULL; AliTRDv0Info *v0(NULL); for(Int_t itrk=0; itrkGetEntriesFast(); itrk++){ track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk); for(Int_t iv(0); ivGetEntriesFast(); iv++){ if(!(v0 = (AliTRDv0Info*)fV0s->At(iv))) continue; if(!v0->HasTrack(track)) continue; ULong_t status = track->GetStatus(); if(!(status&AliESDtrack::kTRDpid)) continue; fhQualityReductions->Fill(v0->GetQuality());//fills integer codes for tracks cut out by track/V0 quality cuts if(!(v0->GetQuality() == 1)) continue; for(Int_t part = 0; part < AliPID::kSPECIES; part++){ fhCutReductions[part]->Fill(v0->GetPID(part,track));//fill in numbers of tracks eliminated by different PID cuts } for(Int_t idecay(0), part(-1); idecay < Int_t(AliTRDv0Info::kNDecays); idecay++){//loop over decay types considered for reference data switch(idecay){ case AliTRDv0Info::kLambda: //protons and pions from Lambda case AliTRDv0Info::kAntiLambda: //antiprotons and pions from Anti-Lambda part = AliPID::kProton; break; case AliTRDv0Info::kK0s: //pions from K0s part = AliPID::kPion; break; case AliTRDv0Info::kGamma: //electrons from conversions part = AliPID::kElectron; break; } //fill histograms with track/V0 quality cuts only fhPsiPair[idecay][0]->Fill(v0->GetV0Momentum(),v0->GetPsiPair());//Angle between daughter momentum plane and plane perpendicular to magnetic field fhInvMass[idecay]->Fill(v0->GetV0Momentum(),v0->GetInvMass(idecay));//Invariant mass fhPointAngle[idecay][0]->Fill(v0->GetV0Momentum(),v0->GetPointingAngle());// = TMath::ACos(esdv0->GetV0CosineOfPointingAngle()); // Cosine of pointing angle fhOpenAngle[idecay][0]->Fill(v0->GetV0Momentum(),v0->GetOpenAngle());// opening angle between daughters fhDCA[idecay][0]->Fill(v0->GetV0Momentum(),v0->GetDCA());// Distance of closest approach of daughter tracks fhV0Chi2ndf[idecay][0]->Fill(v0->GetV0Momentum(),v0->GetChi2ndf(idecay));//Kalman Filter Chi2/NDF fhRadius[idecay][0]->Fill(v0->GetV0Momentum(),v0->GetRadius());//distance of decay/conversion from primary vertex in x-y plane if(v0->HasTrack(track) == -1){ fhTPCdEdx[part][0]->Fill(v0->GetV0Daughter(-1)->P(),v0->GetTPCdEdx(AliTRDv0Info::kNeg));//TPC dE/dx for negative track } else if(v0->HasTrack(track) == 1){ fhTPCdEdx[part][0]->Fill(v0->GetV0Daughter(1)->P(),v0->GetTPCdEdx(AliTRDv0Info::kPos));//TPC dE/dx for positive track } //fill histograms after invariant mass cuts if((v0->GetInvMass(idecay) < v0->GetUpInvMass(idecay,0))&&(v0->GetInvMass(idecay)> v0->GetDownInvMass(idecay))){ fhV0Chi2ndf[idecay][1]->Fill(v0->GetV0Momentum(),v0->GetChi2ndf(idecay)); fhPsiPair[idecay][1]->Fill(v0->GetV0Momentum(),v0->GetPsiPair()); fhPointAngle[idecay][1]->Fill(v0->GetV0Momentum(),v0->GetPointingAngle()); fhOpenAngle[idecay][1]->Fill(v0->GetV0Momentum(),v0->GetOpenAngle()); fhDCA[idecay][1]->Fill(v0->GetV0Momentum(),v0->GetDCA()); fhRadius[idecay][1]->Fill(v0->GetV0Momentum(),v0->GetRadius()); if(v0->HasTrack(track) == -1) fhTPCdEdx[part][1]->Fill(v0->GetV0Daughter(-1)->P(),v0->GetTPCdEdx(AliTRDv0Info::kNeg)); else if(v0->HasTrack(track) == 1) fhTPCdEdx[part][1]->Fill(v0->GetV0Daughter(1)->P(),v0->GetTPCdEdx(AliTRDv0Info::kPos)); } //fill histograms after all reference selection cuts if(v0->GetPID(part,track)==1){ fhV0Chi2ndf[idecay][2]->Fill(v0->GetV0Momentum(),v0->GetChi2ndf(idecay)); fhPsiPair[idecay][2]->Fill(v0->GetV0Momentum(),v0->GetPsiPair()); fhPointAngle[idecay][2]->Fill(v0->GetV0Momentum(),v0->GetPointingAngle()); fhOpenAngle[idecay][2]->Fill(v0->GetV0Momentum(),v0->GetOpenAngle()); fhDCA[idecay][2]->Fill(v0->GetV0Momentum(),v0->GetDCA()); fhRadius[idecay][2]->Fill(v0->GetV0Momentum(),v0->GetRadius()); if(v0->HasTrack(track) == -1) fhTPCdEdx[part][2]->Fill(v0->GetV0Daughter(-1)->P(),v0->GetTPCdEdx(AliTRDv0Info::kNeg)); else if(v0->HasTrack(track) == 1) fhTPCdEdx[part][2]->Fill(v0->GetV0Daughter(1)->P(),v0->GetTPCdEdx(AliTRDv0Info::kPos)); } } } } }