/************************************************************************** * 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-commercial 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. * **************************************************************************/ //------------------------------------------------------------------------------ // Implementation of AliPerformancePtCalib class. It fills histograms with ESD or // TPC track information to study charge/pt spectra. // To analyse the output of AliPerformancePtCalib use AliPerfAnalyzeInvPt class: // Projection of charge/pt vs theta and vs phi resp. histoprams will be fitted // with either polynomial or gaussian fit function to extract minimum position of // charge/pt. // Fit options and theta, phi bins can be set by user. // Attention: use the Set functions of AliPerformancePtCalib when running // AliPerformancePtCalib::Analyse() // The result of the analysis (histograms/graphs) are stored in the folder which is // a data member of AliPerformancePtCalib. // // Author: S.Schuchmann 11/13/2009 // sschuchm@ikf.uni-frankfurt.de //------------------------------------------------------------------------------ /* // after running the performance task, read the file, and get component TFile f("Output.root"); AliPerformancePtCalib * compObj = (AliPerformancePtCalib*)coutput->FindObject("AliPerformancePtCalib"); // set phi and theta bins for fitting and analyse comparison data compObj->SetProjBinsTheta(thetaBins,nThetaBins,minPhi, maxPhi); compObj->SetProjBinsPhi(phiBins,nPhiBins,minTheta,maxTheta); compObj->SetMakeFitOption(kFALSE,exclRange,fitRange); compObj->SetDoRebin(rebin); compObj->Analyse(); //details see functions of this class // the output histograms/graphs will be stored in the folder "folderRes" compObj->GetAnalysisFolder()->ls("*"); // user can save whole comparison object (or only folder with anlysed histograms) // in the seperate output file (e.g.) TFile fout("Analysed_InvPt.root","recreate"); compObj->Write(); // compObj->GetAnalysisFolder()->Write(); fout.Close(); */ #include "TH1F.h" #include "TH2F.h" #include "THnSparse.h" #include "TList.h" #include "TMath.h" #include "TFolder.h" #include "AliESDEvent.h" #include "AliESDtrack.h" #include "AliESDfriendTrack.h" #include "AliESDfriend.h" #include "AliESDtrackCuts.h" #include "AliESDpid.h" #include "AliPerfAnalyzeInvPt.h" #include "AliPerformancePtCalib.h" using namespace std; ClassImp(AliPerformancePtCalib) //________________________________________________________________________ AliPerformancePtCalib::AliPerformancePtCalib(const Char_t* name, const Char_t* title): AliPerformanceObject(name,title), // option parameter for AliPerformancePtCalib::Analyse() fNThetaBins(0), fNPhiBins(0), fMaxPhi(0), fMinPhi(0), fMaxTheta(0), fMinTheta(0), fRange(0), fExclRange(0), fFitGaus(0) , fDoRebin(0), fRebin(0), // option parameter for user defined charge/pt shift fShift(0), fDeltaInvP(0), //options for cuts fOptTPC(0), fESDcuts(0), fPions(0), fEtaAcceptance(0), fCutsRC(0), fCutsMC(0), fList(0), // histograms fHistInvPtPtThetaPhi(0), fHistPtShift0(0), fHistPrimaryVertexPosX(0), fHistPrimaryVertexPosY(0), fHistPrimaryVertexPosZ(0), fHistTrackMultiplicity(0), fHistTrackMultiplicityCuts(0), fHistTPCMomentaPosP(0), fHistTPCMomentaNegP(0), fHistTPCMomentaPosPt(0), fHistTPCMomentaNegPt(0), fHistUserPtShift(0), fHistdedxPions(0), //esd track cuts fESDTrackCuts(0), //pid fESDpid(0), // analysis folder fAnalysisFolder(0) { // Constructor fShift = kFALSE; fDeltaInvP = 0.00; //options for cuts fOptTPC = kTRUE; // read TPC tracks yes/no fESDcuts = kFALSE; fPions = kFALSE; fCutsRC = NULL; fCutsMC = NULL; //esd track cut options fEtaAcceptance = 0.8; // options for function AliPerformancePtCalibMC::Analyse() fFitGaus = kFALSE;// use gaussian function for fitting charge/pt yes/no fNThetaBins = 0; //number of theta bins fNPhiBins = 0; //number of phi bins fMaxPhi = 6.5;// max phi for 2D projection on theta and charge/pt axis fMinPhi = 0.0;// min phi for 2D projection on theta and charge/pt axis fMaxTheta = 3.0;// max theta for 2D projection on phi and charge/pt axis fMinTheta = 0.0;// min theta for 2D projection on phi and charge/pt axis fRange = 0; //fit range around 0 fExclRange =0; //range of rejection of points around 0 fDoRebin = kFALSE; fRebin = 0; for(Int_t i=0; i<100; i++) { fThetaBins[i] = 0.0; fPhiBins[i] = 0.0; } Init(); } //________________________________________________________________________ AliPerformancePtCalib::~AliPerformancePtCalib() { // // destructor // if(fList) delete fList; // histograms if( fHistInvPtPtThetaPhi) delete fHistInvPtPtThetaPhi;fHistInvPtPtThetaPhi=0; if(fHistPtShift0) delete fHistPtShift0;fHistPtShift0=0; if(fHistPrimaryVertexPosX) delete fHistPrimaryVertexPosX;fHistPrimaryVertexPosX=0; if(fHistPrimaryVertexPosY) delete fHistPrimaryVertexPosY;fHistPrimaryVertexPosY=0; if(fHistPrimaryVertexPosZ) delete fHistPrimaryVertexPosZ;fHistPrimaryVertexPosZ=0; if(fHistTrackMultiplicity) delete fHistTrackMultiplicity;fHistTrackMultiplicity=0; if(fHistTrackMultiplicityCuts) delete fHistTrackMultiplicityCuts;fHistTrackMultiplicityCuts=0; if(fHistTPCMomentaPosP) delete fHistTPCMomentaPosP;fHistTPCMomentaPosP=0; if(fHistTPCMomentaNegP) delete fHistTPCMomentaNegP;fHistTPCMomentaNegP=0; if(fHistTPCMomentaPosPt) delete fHistTPCMomentaPosPt;fHistTPCMomentaPosPt=0; if(fHistTPCMomentaNegPt) delete fHistTPCMomentaNegPt ;fHistTPCMomentaNegPt=0; if(fHistUserPtShift) delete fHistUserPtShift;fHistUserPtShift=0; if(fHistdedxPions) delete fHistdedxPions;fHistdedxPions=0; //esd track cuts if(fESDTrackCuts) delete fESDTrackCuts;fESDTrackCuts=0; //pid if(fESDpid) delete fESDpid;fESDpid=0; //analysis folder if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0; } //________________________________________________________________________ void AliPerformancePtCalib::Init() { // Create histograms // Called once fList = new TList(); // init folder fAnalysisFolder = CreateFolder("folderPt_TPC","Analysis Pt Resolution Folder"); // Primary Vertex: fHistPrimaryVertexPosX = new TH1F("fHistPrimaryVertexPosX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5); fList->Add(fHistPrimaryVertexPosX); fHistPrimaryVertexPosY = new TH1F("fHistPrimaryVertexPosY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5); fList->Add(fHistPrimaryVertexPosY); fHistPrimaryVertexPosZ = new TH1F("fHistPrimaryVertexPosZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-2.0,2.0); fList->Add(fHistPrimaryVertexPosZ); // Multiplicity: fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 250, 0, 250); fList->Add(fHistTrackMultiplicity); fHistTrackMultiplicityCuts = new TH1F("fHistTrackMultiplicityCuts", "Multiplicity distribution;Number of tracks after cuts;Events", 250, 0, 250); fList->Add(fHistTrackMultiplicityCuts); // momentum histos //pt shift 0 only needed if shift in 1/pt is applied fHistPtShift0 = new TH1F("fHistPtShift0","1/pt dN/pt vs. pt of ESD track ",800,-40.0,40.0); fList->Add(fHistPtShift0); // THnSparse for 1/pt and pt spectra vs angles const Int_t invPtDims = 4; fMaxPhi = 6.52; fMinPhi = 0.0; fMaxTheta = 3.0; fMinTheta = 0.0; Double_t xminInvPt[invPtDims] = {-4.5,-40.0,fMinTheta,fMinPhi}; Double_t xmaxInvPt[invPtDims] = {4.5,40.0,fMaxTheta,fMaxPhi}; Int_t binsInvPt[invPtDims] = {450,400,150,163}; fHistInvPtPtThetaPhi = new THnSparseF("fHistInvPtPtThetaPhi","1/pt vs pt vs #theta vs #phi ",invPtDims,binsInvPt,xminInvPt,xmaxInvPt); fList->Add(fHistInvPtPtThetaPhi); // momentum test histos fHistTPCMomentaPosP = new TH2F("fHistTPCMomentaPosP","TPC p vs global esd track p pos",300,0.0,15.0,300,0.0,15.0); fList->Add(fHistTPCMomentaPosP); fHistTPCMomentaNegP = new TH2F("fHistTPCMomentaNegP","TPC p vs global esd track p neg",300,0.0,15.0,300,0.0,15.0); fList->Add(fHistTPCMomentaNegP); fHistTPCMomentaPosPt = new TH2F("fHistTPCMomentaPosPt","TPC pt vs global esd track pt pos",300,0.0,15.0,300,0.0,15.0); fList->Add(fHistTPCMomentaPosPt); fHistTPCMomentaNegPt = new TH2F("fHistTPCMomentaNegPt","TPC pt vs global esd track pt neg",300,0.0,15.0,300,0.0,15.0); fList->Add(fHistTPCMomentaNegPt); //user pt shift check fHistUserPtShift = new TH1F("fHistUserPtShift","user defined shift in 1/pt",100,-0.5,1.5); fList->Add(fHistUserPtShift); //pid by dedx fHistdedxPions = new TH2F ("fHistdedxPions","dEdx of pions ident. via kPID vs signed Pt",300,-15.05,15.05,200,0.0,400.0); fList->Add(fHistdedxPions); //pid fESDpid = new AliESDpid(); // esd track cuts fESDTrackCuts =NULL; } //________________________________________________________________________ void AliPerformancePtCalib::SetPtShift(const Double_t shiftVal ) { //set user defined shift in charge/pt if(shiftVal) { fShift=kTRUE; fDeltaInvP = shiftVal; } } //________________________________________________________________________ void AliPerformancePtCalib::Exec(AliMCEvent* const /*mcEvent*/, AliESDEvent *const esdEvent, AliESDfriend * const /*esdFriend*/, const Bool_t /*bUseMC*/, const Bool_t /*bUseESDfriend*/) { //exec: read esd or tpc if(!fESDTrackCuts) { Printf("no esd track cut"); return; } if (!esdEvent) { Printf("ERROR: Event not available"); return; } if (!(esdEvent->GetNumberOfTracks())) return; //vertex info for cut const AliESDVertex *vtx = esdEvent->GetPrimaryVertex(); if (!vtx->GetStatus()) return ; //histo fo user defined shift in charge/pt if(fShift) fHistUserPtShift->Fill(fDeltaInvP); //trakc multiplicity fHistTrackMultiplicity->Fill(esdEvent->GetNumberOfTracks()); // read primary vertex info Double_t tPrimaryVtxPosition[3]; const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexTPC(); tPrimaryVtxPosition[0] = primaryVtx->GetXv(); tPrimaryVtxPosition[1] = primaryVtx->GetYv(); tPrimaryVtxPosition[2] = primaryVtx->GetZv(); fHistPrimaryVertexPosX->Fill(tPrimaryVtxPosition[0]); fHistPrimaryVertexPosY->Fill(tPrimaryVtxPosition[1]); fHistPrimaryVertexPosZ->Fill(tPrimaryVtxPosition[2]); //_fill histos for pt spectra and shift of transverse momentum Int_t count=0; for(Int_t j = 0;jGetNumberOfTracks();j++){// track loop AliESDtrack *esdTrack = esdEvent->GetTrack(j); if(!esdTrack) continue; if(fESDcuts){ if(!fESDTrackCuts->AcceptTrack(esdTrack))continue; } // fill histos if(fOptTPC){ //TPC tracks const AliExternalTrackParam *tpcTrack = esdTrack->GetTPCInnerParam(); if(!tpcTrack) continue; if(fabs(tpcTrack->Eta())>= fEtaAcceptance) continue; Double_t signedPt = tpcTrack->GetSignedPt(); // pid if(fPions){ fESDpid->GetTPCResponse().SetBetheBlochParameters(0.0283086,2.63394e+01,5.04114e-11, 2.12543e+00,4.88663e+00); if( TMath::Abs(fESDpid->NumberOfSigmasTPC(esdTrack,AliPID::kPion)) >1) continue; fHistdedxPions->Fill(signedPt,esdTrack->GetTPCsignal()); } Double_t invPt = 0.0; if(signedPt) { invPt = 1.0/signedPt; fHistPtShift0->Fill(signedPt); if(fShift){Printf("user shift of momentum SET to non zero value!"); invPt += fDeltaInvP; //shift momentum for tests if(invPt) signedPt = 1.0/invPt; else continue; } Double_t theta = tpcTrack->Theta(); Double_t phi = tpcTrack->Phi(); Double_t momAng[4] = {invPt,signedPt,theta,phi}; fHistInvPtPtThetaPhi->Fill(momAng); Double_t pTPC = tpcTrack->GetP(); Double_t pESD = esdTrack->GetP(); Double_t ptESD = esdTrack->GetSignedPt(); if(esdTrack->GetSign()>0){ //compare momenta ESD track and TPC track fHistTPCMomentaPosP->Fill(fabs(pESD),fabs(pTPC)); fHistTPCMomentaPosPt->Fill(fabs(ptESD),fabs(signedPt)); } else{ fHistTPCMomentaNegP->Fill(fabs(pESD),fabs(pTPC)); fHistTPCMomentaNegPt->Fill(fabs(ptESD),fabs(signedPt)); } count++; } else continue; } else{// ESD tracks if(fabs(esdTrack->Eta())> fEtaAcceptance) continue; Double_t invPt = 0.0; Double_t signedPt = esdTrack->GetSignedPt(); if(signedPt){ invPt = 1.0/signedPt; fHistPtShift0->Fill(signedPt); if(fShift){Printf("user shift of momentum SET to non zero value!"); invPt += fDeltaInvP;//shift momentum for tests if(invPt) signedPt = 1.0/invPt; else continue; } Double_t theta = esdTrack->Theta(); Double_t phi = esdTrack->Phi(); Double_t momAng[4] = {invPt,signedPt,theta,phi}; fHistInvPtPtThetaPhi->Fill(momAng); count++; } } } fHistTrackMultiplicityCuts->Fill(count); } //______________________________________________________________________________________________________________________ void AliPerformancePtCalib::Analyse() { // analyse charge/pt spectra in bins of theta and phi. Bins can be set by user THnSparseF *copyTHnSparseTheta = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparseTheta"); if(!copyTHnSparseTheta) return; copyTHnSparseTheta->GetAxis(3)->SetRangeUser(fMinPhi,fMaxPhi); TH2F *histInvPtTheta = (TH2F*)copyTHnSparseTheta->Projection(2,0); THnSparseF *copyTHnSparsePhi = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparsePhi"); if(!copyTHnSparsePhi) return; copyTHnSparsePhi->GetAxis(2)->SetRangeUser(fMinTheta,fMaxTheta); TH2F *histInvPtPhi = (TH2F*)copyTHnSparsePhi->Projection(3,0); AliPerfAnalyzeInvPt *ana = new AliPerfAnalyzeInvPt("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt"); if(!ana) return; TH1::AddDirectory(kFALSE); ana->SetProjBinsTheta(fThetaBins,fNThetaBins); ana->SetProjBinsPhi(fPhiBins,fNPhiBins); ana->SetMakeFitOption(fFitGaus,fExclRange,fRange); if(fDoRebin) ana->SetDoRebin(fRebin); TObjArray *aFolderObj = new TObjArray; if(!aFolderObj) return; ana->StartAnalysis(histInvPtTheta,histInvPtPhi,aFolderObj); // export objects to analysis folder fAnalysisFolder = ExportToFolder(aFolderObj); // delete only TObjArray if(aFolderObj) delete aFolderObj; if(ana) delete ana; } //______________________________________________________________________________________________________________________ TFolder* AliPerformancePtCalib::ExportToFolder(TObjArray * array) { // recreate folder every time and export objects to new one // AliPerformancePtCalib * comp=this; TFolder *folder = comp->GetAnalysisFolder(); TString name, title; TFolder *newFolder = 0; Int_t i = 0; Int_t size = array->GetSize(); if(folder) { // get name and title from old folder name = folder->GetName(); title = folder->GetTitle(); // delete old one delete folder; // create new one newFolder = CreateFolder(name.Data(),title.Data()); newFolder->SetOwner(); // add objects to folder while(i < size) { newFolder->Add(array->At(i)); i++; } } return newFolder; } //______________________________________________________________________________________________________________________ Long64_t AliPerformancePtCalib::Merge(TCollection* const list) { // Merge list of objects (needed by PROOF) if (!list) return 0; if (list->IsEmpty()) return 1; TIterator* iter = list->MakeIterator(); TObject* obj = 0; // collection of generated histograms Int_t count=0; while((obj = iter->Next()) != 0) { AliPerformancePtCalib* entry = dynamic_cast(obj); if (entry == 0) continue; fHistInvPtPtThetaPhi->Add(entry->fHistInvPtPtThetaPhi); fHistPtShift0->Add(entry->fHistPtShift0); fHistPrimaryVertexPosX->Add(entry->fHistPrimaryVertexPosX); fHistPrimaryVertexPosY->Add(entry->fHistPrimaryVertexPosY); fHistPrimaryVertexPosZ->Add(entry->fHistPrimaryVertexPosZ); fHistTrackMultiplicity->Add(entry->fHistTrackMultiplicity); fHistTrackMultiplicityCuts->Add(entry->fHistTrackMultiplicityCuts); fHistTPCMomentaPosP->Add(entry->fHistTPCMomentaPosP); fHistTPCMomentaNegP->Add(entry->fHistTPCMomentaNegP); fHistTPCMomentaPosPt->Add(entry->fHistTPCMomentaPosPt); fHistTPCMomentaNegPt->Add(entry->fHistTPCMomentaNegPt); fHistdedxPions->Add(entry->fHistdedxPions); count++; } return count; } //______________________________________________________________________________________________________________________ TFolder* AliPerformancePtCalib::CreateFolder(TString name,TString title) { // create folder for analysed histograms // TFolder *folder = 0; folder = new TFolder(name.Data(),title.Data()); return folder; } //______________________________________________________________________________________________________________________ void AliPerformancePtCalib::SetProjBinsPhi(const Double_t *phiBinArray,const Int_t nphBins, const Double_t minTheta, const Double_t maxTheta){ // set phi bins for Analyse() //set phi bins as array and set number of this array which is equal to number of bins analysed //the last analysed bin will always be the projection from first to last bin in the array if(nphBins){ fNPhiBins = nphBins; for(Int_t k = 0;k