From b168fbc0c6c6d8308bef9070e911c99ed23ef18f Mon Sep 17 00:00:00 2001 From: marian Date: Mon, 6 Aug 2007 10:28:11 +0000 Subject: [PATCH] Up to date of the components (Marian, Lars, Alexander) --- TPC/TPCcalib/AliTPCSelectorTracks.cxx | 72 +- TPC/TPCcalib/AliTPCcalibTracks.cxx | 2541 +++++++++++++++++++------ TPC/TPCcalib/AliTPCcalibTracks.h | 231 ++- TPC/TPCcalib/AliTPCcalibV0.cxx | 308 ++- TPC/TPCcalib/AliTPCcalibV0.h | 4 +- TPC/TPCcalib/make_TPCcalib_par.sh | 1 + TPC/TPCcalib/tpcSelectorTracks.C | 2 +- 7 files changed, 2528 insertions(+), 631 deletions(-) diff --git a/TPC/TPCcalib/AliTPCSelectorTracks.cxx b/TPC/TPCcalib/AliTPCSelectorTracks.cxx index cbfc229154e..632baeda9b0 100644 --- a/TPC/TPCcalib/AliTPCSelectorTracks.cxx +++ b/TPC/TPCcalib/AliTPCSelectorTracks.cxx @@ -30,8 +30,8 @@ #include "TCint.h" #include "TH1I.h" // -#include "AliMagF.h" #include "AliTracker.h" +#include "AliMagF.h" // #include "AliESD.h" #include "AliESDtrack.h" @@ -42,6 +42,7 @@ #include "AliClusterMap.h" // #include "AliTPCcalibTracks.h" +//#include "AliTPCcalibTracks.cxx" #include "AliTPCSelectorTracks.h" @@ -82,10 +83,9 @@ void AliTPCSelectorTracks::SlaveBegin(TTree * tree) fOutput->AddLast(fNtracks); fOutput->AddLast(fNtracksFriend); fOutput->AddLast(fNClusters); - fCalibTracks = new AliTPCcalibTracks; - // - fCalibTracks->ProofSlaveBegin(fOutput); - + + fCalibTracks = new AliTPCcalibTracks("calibTracks", "calibTracks"); + fOutput->AddLast(fCalibTracks); } @@ -127,30 +127,30 @@ Bool_t AliTPCSelectorTracks::Process(Long64_t entry) try { fChain->GetTree()->GetEntry(entry); } catch (std::bad_alloc) { - printf("Pica vyjebana pojebany skurveny kokot piciak\n"); + Info("Procces","Pica vyjebana pojebany skurveny kokot piciak\n"); fESD =0; fESDfriend = 0; return 0; } // - //Info("Procces","0"); + // Info("Procces","0"); if (!fESD) { - fESD =0; - fESDfriend=0; + fESD = 0; + fESDfriend = 0; //CleanESD(); return kFALSE; } Int_t ntracks = fESD->GetNumberOfTracks(); fNtracks->Fill(ntracks); - //Info("Procces",Form("1-Ntracks=%d",ntracks)); + Info("Procces", Form("entry\t%d: Ntracks = %d",entry, ntracks)); - if (!fESDfriend || fESDfriend->GetNumberOfTracks()!=ntracks) { - fESD =0; - fESDfriend=0; + if (!fESDfriend || fESDfriend->GetNumberOfTracks() != ntracks) { + fESD = 0; + fESDfriend = 0; // CleanESD(); if (fESDfriend) fNtracksFriend->Fill(fESDfriend->GetNumberOfTracks()); - //Info("Procces","2- PROBLEM"); + Info("Procces","2: PROBLEM"); return kFALSE; } fESD->SetESDfriend(fESDfriend); @@ -160,10 +160,11 @@ Bool_t AliTPCSelectorTracks::Process(Long64_t entry) AliTPCseed *seed; AliTPCclusterMI cl; - for (Int_t tr=0; tr < ntracks; tr++){ - AliESDtrack * esdTrack = (AliESDtrack*) fESD->GetTrack(tr); - AliESDfriendTrack *friendtrack = (AliESDfriendTrack*) fESD->GetTrack(tr)->GetFriendTrack(); + for (Int_t tr = 0; tr < ntracks; tr++){ + AliESDtrack *esdTrack = (AliESDtrack*) fESD->GetTrack(tr); + AliESDfriendTrack *friendtrack = (AliESDfriendTrack*) fESD->GetTrack(tr)->GetFriendTrack(); seed = (AliTPCseed*)(friendtrack->GetCalibObject(0)); + Info("Process",Form("Proccessing track%d\n",tr)); if (seed) { if (!fCalibTracks->AcceptTrack(seed)) continue; // FillHistoCluster(seed); @@ -174,7 +175,7 @@ Bool_t AliTPCSelectorTracks::Process(Long64_t entry) // } } - CleanESD(); + CleanESD(); return kTRUE; } @@ -194,11 +195,26 @@ void AliTPCSelectorTracks::Terminate() // The Terminate() function is the last function to be called during // a query. It always runs on the client, it can be used to present // the results graphically or save the results to file. - + + printf ("Terminate... \n"); if (!fOutput) return; TFile file("Output.root","recreate"); + printf("fOutput contains the following: \n"); + fOutput->Print(); + printf("Trying to write the file 'Output.root'... \n"); fOutput->Write(); + file.Close(); + +// fCalibTracks->MakePlots(1); +// fCalibTracks->FitResolutionNew(1); +// fCalibTracks->MakeReport(1); + + + + } + + void AliTPCSelectorTracks::Init(TTree *tree) { // The Init() function is called when the selector needs to initialize @@ -216,15 +232,15 @@ void AliTPCSelectorTracks::Init(TTree *tree) tree->SetBranchStatus("*",1); // fChain->SetMakeClass(1); fChain->SetBranchAddress("ESD",&fESD); - //Info("Init","Enter"); + Info("Init","Enter"); Bool_t isOK=kFALSE; if (fChain->GetBranch("ESDfriend")) { fChain->SetBranchAddress("ESDfriend",&fESDfriend); - //Info("Init","V0-ESDfriend."); + Info("Init","V0-ESDfriend."); isOK=kTRUE; } if (fChain->GetBranch("ESDfriend.")){ - //Info("Init","V1-ESDfriend."); + Info("Init","V1-ESDfriend."); fChain->SetBranchAddress("ESDfriend.",&fESDfriend); isOK=kTRUE; } @@ -234,21 +250,21 @@ void AliTPCSelectorTracks::Init(TTree *tree) // Try to solve problem // - //Info("Init","Problem"); + Info("Init","Problem"); if (tree->GetBranch("ESD")){ - //Info("InitTree",tree->GetBranch("ESD")->GetFile()->GetName()); + Info("InitTree",tree->GetBranch("ESD")->GetFile()->GetName()); char fname[1000]; sprintf(fname,"%s/AliESDfriends.root",gSystem->DirName(tree->GetBranch("ESD")->GetFile()->GetName())); - //Info("InitFile",fname); + Info("InitFile",fname); if (tree->AddFriend("esdFriendTree",fname)){ - //Info("InitFileOK",fname); + Info("InitFileOK",fname); if (fChain->GetBranch("ESDfriend")) { fChain->SetBranchAddress("ESDfriend",&fESDfriend); - //Info("Init","V0-ESDfriend."); + Info("Init","V0-ESDfriend."); isOK=kTRUE; } if (fChain->GetBranch("ESDfriend.")){ - //Info("Init","V1-ESDfriend."); + Info("Init","V1-ESDfriend."); fChain->SetBranchAddress("ESDfriend.",&fESDfriend); isOK=kTRUE; } diff --git a/TPC/TPCcalib/AliTPCcalibTracks.cxx b/TPC/TPCcalib/AliTPCcalibTracks.cxx index 79c7c31802d..4d322882b6e 100644 --- a/TPC/TPCcalib/AliTPCcalibTracks.cxx +++ b/TPC/TPCcalib/AliTPCcalibTracks.cxx @@ -1,3 +1,40 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + + +/////////////////////////////////////////////////////////////////////////////// +// // +// Class to analyse tracks for calibration // +// to be used as a component in AliTPCSelectorTracks // +// In the constructor you have to specify name and title // +// to get the Object out of a file. // +// The parameter 'clusterParam', a AliTPCClusterParam object // +// (needed for TPC cluster error and shape parameterization) // +// Normally you get this object out of the file 'TPCClusterParam.root' // +// In the parameter 'cuts' the cuts are specified, that decide // +// weather a track will be accepted for calibration or not. // +// // +// // +/////////////////////////////////////////////////////////////////////////////// + +// +// ROOT includes +// +#include +#include +using namespace std; #include #include @@ -10,9 +47,14 @@ #include "TMatrixD.h" #include "TTreeStream.h" #include "TF1.h" +#include +#include +#include "TPostScript.h" +#include "TCint.h" - - +// +// AliROOT includes +// #include "AliMagF.h" #include "AliTracker.h" #include "AliESD.h" @@ -23,53 +65,291 @@ #include "AliTPCclusterMI.h" #include "AliTPCROC.h" - #include "AliTPCParamSR.h" -#include "AliTPCClusterParam.h" #include "AliTrackPointArray.h" -#include "TCint.h" #include "AliTPCcalibTracks.h" +#include "AliTPCClusterParam.h" + ClassImp(AliTPCcalibTracks) AliTPCParam param; +AliTPCcalibTracks::AliTPCcalibTracks():TNamed() { + // + // AliTPCcalibTracks default constructor + // + cout << "AliTPCcalibTracks' default constructor called" << endl; + fClusterParam = 0; + fArrayAmpRow = 0; + fArrayAmp = 0; + fArrayQDY = 0; + fArrayQDZ = 0; + fArrayQRMSY = 0; + fArrayQRMSZ = 0; + fDeltaY = 0; + fDeltaZ = 0; + fResolY = 0; + fResolZ = 0; + fRMSY = 0; + fRMSZ = 0; + fHclus = 0; + fROC = 0; + fCuts = 0; +} + -AliTPCcalibTracks::AliTPCcalibTracks() : - TNamed(), - fHclus(0), - fFileNo(0) +AliTPCcalibTracks::AliTPCcalibTracks(AliTPCcalibTracks* ct){ + // + // AliTPCcalibTracks copy constructor + // + + cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl; + + Int_t length = ct->fArrayAmpRow->GetEntriesFast(); + fArrayAmpRow = new TObjArray(length); + fArrayAmp = new TObjArray(length); + for (Int_t i = 0; i < length; i++) { + fArrayAmpRow->AddAt( (TProfile*)ct->fArrayAmpRow->At(i)->Clone(), i); + fArrayAmp->AddAt( ((TProfile*)ct->fArrayAmp->At(i)->Clone()), i); + } + + length = ct->fArrayQDY->GetEntriesFast(); + fArrayQDY= new TObjArray(length); + fArrayQDZ= new TObjArray(length); + fArrayQRMSY= new TObjArray(length); + fArrayQRMSZ= new TObjArray(length); + for (Int_t i = 0; i < length; i++) { + fArrayQDY->AddAt( ((TH1F*)ct->fArrayQDY->At(i)->Clone()), i); + fArrayQDZ->AddAt( ((TH1F*)ct->fArrayQDZ->At(i)->Clone()), i); + fArrayQRMSY->AddAt( ((TH1F*)ct->fArrayQRMSY->At(i)->Clone()), i); + fArrayQRMSZ->AddAt( ((TH1F*)ct->fArrayQRMSZ->At(i)->Clone()), i); + } + + length = ct->fResolY->GetEntriesFast(); + fResolY = new TObjArray(length); + fResolZ = new TObjArray(length); + fRMSY = new TObjArray(length); + fRMSZ = new TObjArray(length); + for (Int_t i = 0; i < length; i++) { + fResolY->AddAt( ((TH1F*)ct->fResolY->At(i)->Clone()), i); + fResolZ->AddAt( ((TH1F*)ct->fResolZ->At(i)->Clone()), i); + fRMSY->AddAt( ((TH1F*)ct->fRMSY->At(i)->Clone()), i); + fRMSZ->AddAt( ((TH1F*)ct->fRMSZ->At(i)->Clone()), i); + } + + fDeltaY = (TH1F*)ct->fDeltaY->Clone(); + fDeltaZ = (TH1F*)ct->fDeltaZ->Clone(); + + fHclus = (TH1I*)ct->fHclus->Clone(); + + fCuts = new AliTPCcalibTracksCuts(ct->fCuts->GetMinClusters(), ct->fCuts->GetMinRatio(), + ct->fCuts->GetMax1pt(), ct->fCuts->GetEdgeYXCutNoise(), ct->fCuts->GetEdgeThetaCutNoise()); +} + + +AliTPCcalibTracks::~AliTPCcalibTracks() { + // + // AliTPCcalibTracks destructor + // + cout << "AliTPCcalibTracks' destuctor called." << endl; + Int_t length = 0; + if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast(); + for (Int_t i = 0; i < length; i++){ + delete fArrayAmpRow->At(i); + delete fArrayAmp->At(i); + } + delete fArrayAmpRow; + delete fArrayAmp; + + delete fDeltaY; + delete fDeltaZ; + + if (fResolY) length = fResolY->GetEntriesFast(); + for (Int_t i = 0; i < length; i++){ + delete fResolY->At(i); + delete fResolZ->At(i); + delete fRMSY->At(i); + delete fRMSZ->At(i); + } + delete fResolY; + delete fResolZ; + delete fRMSY; + delete fRMSZ; + + if (fArrayQDY) length = fArrayQDY->GetEntriesFast(); + for (Int_t i = 0; i < length; i++){ + delete fArrayQDY->At(i); + delete fArrayQDZ->At(i); + delete fArrayQRMSY->At(i); + delete fArrayQRMSZ->At(i); + } + delete fArrayQDY; + delete fArrayQDZ; + delete fArrayQRMSY; + delete fArrayQRMSZ; + + delete fHclus; +// delete fDebugStream; +} + + +AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts) : + TNamed(name, title), + fHclus(0) + // + // AliTPCcalibTracks constructor + // specify 'name' and 'title' of your object + // specify 'clusterParam', (needed for TPC cluster error and shape parameterization) + // In the parameter 'cuts' the cuts are specified, that decide + // weather a track will be accepted for calibration or not. + // + // All histograms are instatiated in this constructor. + // { G__SetCatchException(0); param.Update(); - TFile fparam("/u/miranov/TPCClusterParam.root"); - fClusterParam = (AliTPCClusterParam *) fparam.Get("Param"); + + fClusterParam = clusterParam; if (fClusterParam){ - //fClusterParam->SetInstance(fClusterParam); - }else{ + fClusterParam->SetInstance(fClusterParam); + } + else { printf("Cluster Param not found\n"); } - fDebugStream = new TTreeSRedirector("TPCSelectorDebug.root"); - } - - -void AliTPCcalibTracks::ProcessTrack(AliTPCseed * seed){ + fCuts = cuts; + fDebugStream = new TTreeSRedirector("TPCSelectorDebug.root"); // needs investigation !!!!! + + TH1::AddDirectory(kFALSE); + + char chname[1000]; + TProfile * prof1=0; + TH1F * his1 =0; + fHclus = new TH1I("hclus","Number of clusters",100,0,200); // valgrind 3 + + // Amplitude - sector - row histograms + fArrayAmpRow = new TObjArray(72); + fArrayAmp = new TObjArray(72); + for (Int_t i = 0; i < 36; i++){ + sprintf(chname,"Amp_row_Sector%d",i); + prof1 = new TProfile(chname,chname,63,0,64); // valgrind 3 193,536 bytes in 354 blocks are still reachable + prof1->SetXTitle("Pad row"); + prof1->SetYTitle("Mean Max amplitude"); + fArrayAmpRow->AddAt(prof1,i); + sprintf(chname,"Amp_row_Sector%d",i+36); + prof1 = new TProfile(chname,chname,96,0,97); // valgrind 3 3,912 bytes in 6 blocks are possibly lost + prof1->SetXTitle("Pad row"); + prof1->SetYTitle("Mean Max amplitude"); + fArrayAmpRow->AddAt(prof1,i+36); + + // amplitude + sprintf(chname,"Amp_Sector%d",i); + his1 = new TH1F(chname,chname,250,0,500); // valgrind + his1->SetXTitle("Max Amplitude (ADC)"); + fArrayAmp->AddAt(his1,i); + sprintf(chname,"Amp_Sector%d",i+36); + his1 = new TH1F(chname,chname,200,0,600); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable + his1->SetXTitle("Max Amplitude (ADC)"); + fArrayAmp->AddAt(his1,i+36); + } + + TH1::AddDirectory(kFALSE); + + fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1); + fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1); + + fResolY = new TObjArray(3); + fResolZ = new TObjArray(3); + fRMSY = new TObjArray(3); + fRMSZ = new TObjArray(3); + TH3F * his3D; + // + his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1); + fResolY->AddAt(his3D,0); + his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1); + fResolY->AddAt(his3D,1); + his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1); + fResolY->AddAt(his3D,2); + // + his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1); + fResolZ->AddAt(his3D,0); + his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1); + fResolZ->AddAt(his3D,1); + his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1); + fResolZ->AddAt(his3D,2); + // + his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8); + fRMSY->AddAt(his3D,0); + his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8); + fRMSY->AddAt(his3D,1); + his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8); + fRMSY->AddAt(his3D,2); + // + his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8); + fRMSZ->AddAt(his3D,0); + his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8); + fRMSZ->AddAt(his3D,1); + his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8); + fRMSZ->AddAt(his3D,2); + // + + TH1::AddDirectory(kFALSE); + + fArrayQDY = new TObjArray(300); + fArrayQDZ = new TObjArray(300); + fArrayQRMSY = new TObjArray(300); + fArrayQRMSZ = new TObjArray(300); + for (Int_t iq = 0; iq <= 10; iq++){ + for (Int_t ipad = 0; ipad < 3; ipad++){ + Int_t bin = GetBin(iq, ipad); + Float_t qmean = GetQ(bin); + char name[200]; + sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1); + fArrayQDY->AddAt(his3D, bin); + sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1); + fArrayQDZ->AddAt(his3D, bin); + sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1); + fArrayQRMSY->AddAt(his3D, bin); + sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1); + fArrayQRMSZ->AddAt(his3D, bin); + } + } +} + + +void AliTPCcalibTracks::Process(AliTPCseed *track, AliESDtrack *esd){ + // + // To be called in the selector + // first AcceptTrack is evaluated, then calls all the following analyse functions: + // FillResolutionHistoLocal(track) + // AlignUpDown(track, esd) + // + if (AcceptTrack(track)) { + FillResolutionHistoLocal(track); + AlignUpDown(track, esd); + } } -Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){ + +Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){ // // calculate bins for given q and pad type // used in TObjArray // - Int_t res = TMath::Max(TMath::Nint((TMath::Sqrt(q)-3.)),0); - res*=3; - res+=pad; + Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 ); + res *= 3; + res += pad; return res; } -Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){ + +Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){ // // calculate bins for given iq and pad type // used in TObjArray @@ -77,187 +357,55 @@ Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){ return iq*3+pad;; } -Float_t AliTPCcalibTracks::GetQ(Int_t bin){ - Int_t bin0 = bin/3; - bin0+=3; - return bin0*bin0; -} - - - - - - - -void AliTPCcalibTracks::ProofSlaveBegin(TList * output) -{ - // Called on PROOF - fill output list - //fChain = tree; - //Init(tree); - - char chname[1000]; - TProfile * prof1=0; - TH1F * his1 =0; - fHclus = new TH1I("hclus","Number of clusters",100,0,200); - output->AddLast(fHclus); - - - // - // Amplitude - sector -row histograms - // - fArrayAmpRow = new TObjArray(72); - fArrayAmp = new TObjArray(72); - for (Int_t i=0; i<36; i++){ - sprintf(chname,"Amp_row_Sector%d",i); - prof1 = new TProfile(chname,chname,63,0,64); - prof1->SetXTitle("Pad row"); - prof1->SetYTitle("Mean Max amplitude"); - fArrayAmpRow->AddAt(prof1,i); - output->AddLast(prof1); - sprintf(chname,"Amp_row_Sector%d",i+36); - prof1 = new TProfile(chname,chname,96,0,97); - prof1->SetXTitle("Pad row"); - prof1->SetYTitle("Mean Max amplitude"); - fArrayAmpRow->AddAt(prof1,i+36); - output->AddLast(prof1); - // - // amplitude - sprintf(chname,"Amp_Sector%d",i); - his1 = new TH1F(chname,chname,250,0,500); - his1->SetXTitle("Max Amplitude (ADC)"); - fArrayAmp->AddAt(his1,i); - output->AddLast(his1); - sprintf(chname,"Amp_Sector%d",i+36); - his1 = new TH1F(chname,chname,200,0,600); - his1->SetXTitle("Max Amplitude (ADC)"); - fArrayAmp->AddAt(his1,i+36); - output->AddLast(his1); - // - } - - fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1); - fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1); - output->AddLast(fDeltaY); - output->AddLast(fDeltaZ); - fResolY = new TObjArray(3); - fResolZ = new TObjArray(3); - fRMSY = new TObjArray(3); - fRMSZ = new TObjArray(3); - TH3F * his3D; - // - his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1); - fResolY->AddAt(his3D,0); - output->AddLast(his3D); - his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1); - fResolY->AddAt(his3D,1); - output->AddLast(his3D); - his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1); - fResolY->AddAt(his3D,2); - output->AddLast(his3D); - // - his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1); - fResolZ->AddAt(his3D,0); - output->AddLast(his3D); - his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1); - fResolZ->AddAt(his3D,1); - output->AddLast(his3D); - his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1); - fResolZ->AddAt(his3D,2); - output->AddLast(his3D); - // - his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8); - fRMSY->AddAt(his3D,0); - output->AddLast(his3D); - his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8); - fRMSY->AddAt(his3D,1); - output->AddLast(his3D); - his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8); - fRMSY->AddAt(his3D,2); - output->AddLast(his3D); - // - his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8); - fRMSZ->AddAt(his3D,0); - output->AddLast(his3D); - his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8); - fRMSZ->AddAt(his3D,1); - output->AddLast(his3D); - his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8); - fRMSZ->AddAt(his3D,2); - output->AddLast(his3D); - // - fArrayQDY = new TObjArray(300); - fArrayQDZ = new TObjArray(300); - fArrayQRMSY = new TObjArray(300); - fArrayQRMSZ = new TObjArray(300); - for (Int_t iq=0; iq<10; iq++){ - for (Int_t ipad=0; ipad<3; ipad++){ - Int_t bin = GetBin(iq,ipad); - Float_t qmean = GetQ(bin); - char name[200]; - sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean); - his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1); - fArrayQDY->AddAt(his3D,bin); - output->AddLast(his3D); - sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean); - his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1); - fArrayQDZ->AddAt(his3D,bin); - output->AddLast(his3D); - // - sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean); - his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1); - fArrayQRMSY->AddAt(his3D,bin); - output->AddLast(his3D); - sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean); - his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1); - fArrayQRMSZ->AddAt(his3D,bin); - output->AddLast(his3D); - } - } +Float_t AliTPCcalibTracks::GetQ(Int_t bin){ + // + // (bin / 3 + 3)^2 + // + Int_t bin0 = bin / 3; + bin0 += 3; + return bin0 * bin0; } - - -Float_t AliTPCcalibTracks::TPCBetheBloch(Float_t bg) -{ - // - // Bethe-Bloch energy loss formula - // - const Double_t kp1=0.76176e-1; - const Double_t kp2=10.632; - const Double_t kp3=0.13279e-4; - const Double_t kp4=1.8631; - const Double_t kp5=1.9479; - Double_t dbg = (Double_t) bg; - Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg); - Double_t aa = TMath::Power(beta,kp4); - Double_t bb = TMath::Power(1./dbg,kp5); - bb=TMath::Log(kp3+bb); - return ((Float_t)((kp2-aa-bb)*kp1/aa)); +Float_t AliTPCcalibTracks::TPCBetheBloch(Float_t bg){ + // + // Bethe-Bloch energy loss formula + // + const Double_t kp1=0.76176e-1; + const Double_t kp2=10.632; + const Double_t kp3=0.13279e-4; + const Double_t kp4=1.8631; + const Double_t kp5=1.9479; + Double_t dbg = (Double_t) bg; + Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg); + Double_t aa = TMath::Power(beta,kp4); + Double_t bb = TMath::Power(1./dbg,kp5); + bb=TMath::Log(kp3+bb); + return ((Float_t)((kp2-aa-bb)*kp1/aa)); } Bool_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){ // + // Function, that decides wheather a given track is accepted for + // the analysis or not. + // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts' // - // - const Int_t kMinClusters = 20; - const Float_t kMinRatio = 0.4; - const Float_t kMax1pt = 0.5; - const Float_t kEdgeYXCutNoise = 0.13; - const Float_t kEdgeThetaCutNoise = 0.018; + const Int_t kMinClusters = fCuts->GetMinClusters(); + const Float_t kMinRatio = fCuts->GetMinRatio(); + const Float_t kMax1pt = fCuts->GetMax1pt(); + const Float_t kEdgeYXCutNoise = fCuts->GetEdgeYXCutNoise(); + const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise(); // // edge induced noise tracks - NEXT RELEASE will be removed during tracking - if (TMath::Abs(track->GetY()/track->GetX())> kEdgeYXCutNoise) - if (TMath::Abs(track->GetTgl())GetNumberOfClusters()GetNumberOfClusters()/(track->GetNFoundable()+1.); - if (ratioGetY() / track->GetX()) > kEdgeYXCutNoise ) + if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return kFALSE; + if (track->GetNumberOfClusters() < kMinClusters) return kFALSE; + Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.); + if (ratio < kMinRatio) return kFALSE; Float_t mpt = track->Get1Pt(); - if (TMath::Abs(mpt)>kMax1pt) return kFALSE; + if (TMath::Abs(mpt) > kMax1pt) return kFALSE; //if (TMath::Abs(track->GetZ())>240.) return kFALSE; //if (TMath::Abs(track->GetZ())<10.) return kFALSE; //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE; @@ -265,171 +413,225 @@ Bool_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){ return kTRUE; } + void AliTPCcalibTracks::FillHistoCluster(AliTPCseed * track){ - // - // - // + // + // fill fArrayAmpRow + // 72 TProfiles, one for each ROC with amplitudes vs. row + // Is this function used somewhere??? + // const Int_t kFirstLargePad = 127; const Float_t kLargePadSize = 1.5; - for (Int_t irow=0; irow<159; irow++){ + for (Int_t irow = 0; irow < 159; irow++){ AliTPCclusterMI * cluster = track->GetClusterPointer(irow); if (!cluster) continue; Int_t sector = cluster->GetDetector(); - if (cluster->GetQ()<=0) continue; + if (cluster->GetQ() <= 0) continue; Float_t max = cluster->GetMax(); - printf ("irow, kFirstLargePad = %d, %d \n",irow,kFirstLargePad); + printf ("irow, kFirstLargePad = %d, %d \n", irow, kFirstLargePad); if ( irow >= kFirstLargePad) { max /= kLargePadSize; } - TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector); + TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector); profAmpRow->Fill(cluster->GetRow(), max); } } + void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ - // - // fill resolution histograms - localy - trcklet in the neighborhood - // - const Int_t kDelta = 10; // delta rows to fit - const Float_t kMinRatio = 0.75; // minimal ratio - const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal - const Float_t kErrorFraction = 0.5; // use only clusters with small intrpolation error - for error param - const Int_t kFirstLargePad = 127; - const Float_t kLargePadSize = 1.5; - static TLinearFitter fitterY2(3,"pol2"); - static TLinearFitter fitterZ2(3,"pol2"); + // + // fill resolution histograms - localy - tracklet in the neighborhood + // write debug information to 'TPCSelectorDebug.root' + // + // + // loop over all padrows along the track + // fit tracklets (length: 13 rows) calculate mean chi^2 for this track-fit in Y and Z direction + // + // loop again over all padrows along the track + // fit tracklet (clusters in given padrow +- kDelta padrows) + // with polynom of 2nd order and two polynoms of 1st order + // take both polynoms of 1st order, calculate difference of their parameters + // add covariance matrixes and calculate chi2 of this difference + // if this chi2 is bigger than a given threshold, assume that the current cluster is + // a kink an goto next padrow + // if not: + // fill fArrayAmpRow, array with amplitudes vs. row for given sector + // fill fArrayAmp, array with amplitude histograms for give sector + // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY + // + // write debug information to 'TPCSelectorDebug.root' + // + + + const Int_t kDelta = 10; // delta rows to fit + const Float_t kMinRatio = 0.75; // minimal ratio + const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal + const Float_t kErrorFraction = 0.5; // use only clusters with small interpolation error - for error param + const Int_t kFirstLargePad = 127; // medium pads -> long pads + const Float_t kLargePadSize = 1.5; // factor between medium and long pads' area static TLinearFitter fitterY0(2,"pol1"); static TLinearFitter fitterZ0(2,"pol1"); static TLinearFitter fitterY1(2,"pol1"); - static TLinearFitter fitterZ1(2,"pol1"); - TVectorD paramY0(2); - TVectorD paramY1(2); - TVectorD paramY2(3); - TVectorD paramZ0(2); - TVectorD paramZ1(2); - TVectorD paramZ2(3); - TMatrixD matrixY0(2,2); - TMatrixD matrixZ0(2,2); - TMatrixD matrixY1(2,2); - TMatrixD matrixZ1(2,2); - // + static TLinearFitter fitterZ1(2,"pol1"); // valgrind 3 20,484 bytes in 435 blocks are indirectly lost + static TLinearFitter fitterY2(3,"pol2"); + static TLinearFitter fitterZ2(3,"pol2"); + TVectorD paramY0(2); + TVectorD paramZ0(2); + TVectorD paramY1(2); + TVectorD paramZ1(2); + TVectorD paramY2(3); + TVectorD paramZ2(3); + TMatrixD matrixY0(2,2); + TMatrixD matrixZ0(2,2); + TMatrixD matrixY1(2,2); + TMatrixD matrixZ1(2,2); + // estimate mean error - // - Int_t nTrackletsAll = 0; - Int_t nClusters = 0; - Float_t csigmaY = 0; - Float_t csigmaZ = 0; + Int_t nTrackletsAll = 0; + Int_t nClusters = 0; + Float_t csigmaY = 0; + Float_t csigmaZ = 0; Int_t sectorG = -1; - for (Int_t irow=0; irow<159; irow++){ + + fHclus->Fill(track->GetNumberOfClusters()); + + for (Int_t irow = 0; irow < 159; irow++){ + // loop over all rows along the track + // fit tracklets (length: 13 rows) with pol2 in Y and Z direction + // calculate mean chi^2 for this track-fit in Y and Z direction AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); - if (!cluster0) continue; + if (!cluster0) continue; // no cluster found Int_t sector = cluster0->GetDetector(); - if (sector!=sectorG){ - nClusters=0; + if (sector != sectorG){ + // track leaves sector before it crossed enough rows to fit / initialization + nClusters = 0; fitterY2.ClearPoints(); fitterZ2.ClearPoints(); - sectorG=sector; - }else{ + sectorG = sector; + } + else { nClusters++; Double_t x = cluster0->GetX(); - fitterY2.AddPoint(&x,cluster0->GetY(),1); - fitterZ2.AddPoint(&x,cluster0->GetZ(),1); + fitterY2.AddPoint(&x, cluster0->GetY(), 1); + fitterZ2.AddPoint(&x, cluster0->GetZ(), 1); // - if (nClusters>=kDelta+3){ + if ( nClusters >= kDelta + 3 ){ + // if more than 13 (kDelta+3) rows / clusters were added to the fitters + // fit the tracklet, increase trackletCounter fitterY2.Eval(); fitterZ2.Eval(); nTrackletsAll++; - csigmaY+=fitterY2.GetChisquare()/(nClusters-3.); - csigmaZ+=fitterZ2.GetChisquare()/(nClusters-3.); - nClusters=-1; + csigmaY += fitterY2.GetChisquare() / (nClusters - 3.); + csigmaZ += fitterZ2.GetChisquare() / (nClusters - 3.); + nClusters = -1; fitterY2.ClearPoints(); fitterZ2.ClearPoints(); } } - } - csigmaY = TMath::Sqrt(csigmaY/nTrackletsAll); - csigmaZ = TMath::Sqrt(csigmaZ/nTrackletsAll); + } // for (Int_t irow = 0; irow < 159; irow++) + // mean chi^2 for all tracklet fits in Y and in Z direction + csigmaY = TMath::Sqrt(csigmaY / nTrackletsAll); + csigmaZ = TMath::Sqrt(csigmaZ / nTrackletsAll); + // // // - for (Int_t irow=0; irow<159; irow++){ + for (Int_t irow = 0; irow < 159; irow++){ + // loop again over all rows along the track + // do analysis + // Int_t nclFound = 0; - Int_t nclFoundable = 0; +// Int_t nclFoundable = 0; AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); if (!cluster0) continue; Int_t sector = cluster0->GetDetector(); Float_t xref = cluster0->GetX(); - // + + + // what is the following loop good for????? +/* // check the neighborhood occupancy - (Delta ray - noise removal) - // - for (Int_t idelta= -kDelta; idelta<=kDelta; idelta++){ - if (idelta==0) continue; - if (idelta+irow<0) continue; - if (idelta+irow>159) continue; + for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ + // loop over irow +- kDelta rows (neighboured rows) + // increase nclFoundable and nclFound + // nclFoundable == nclFound !!!!!!!!! + if (idelta == 0) continue; // check neighbourhood rows, not the row itself + if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC AliTPCclusterMI * clusterD = track->GetClusterPointer(irow); - if ( clusterD && clusterD->GetDetector()!= sector) continue; - if (clusterD->GetType()<0) continue; - nclFoundable++; - if (clusterD) nclFound++; - } - if (nclFoundGetDetector() != sector) continue; // track leaves ROC + if (clusterD->GetType() < 0) continue; + nclFoundable++; // ??????????????????????????????????? + nclFound++; // nclFoundable == nclFound !!!!!!!!! + } // neighbourhood-loop + + if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters found in neighbourhood + if ( Float_t(nclFound) / Float_t(nclFoundable) < kMinRatio ) continue; // if ratio between foundable and found clusters is too bad + +*/ // Make Fit - // fitterY2.ClearPoints(); fitterZ2.ClearPoints(); fitterY0.ClearPoints(); fitterZ0.ClearPoints(); fitterY1.ClearPoints(); fitterZ1.ClearPoints(); - // - nclFound=0; - Int_t ncl0=0; - Int_t ncl1=0; - for (Int_t idelta=-kDelta; idelta<=kDelta; idelta++){ - if (idelta==0) continue; - if (idelta+irow<0) continue; - if (idelta+irow>159) continue; - AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta); + nclFound = 0; + Int_t ncl0 = 0; + Int_t ncl1 = 0; + + // fit tracklet (clusters in given padrow +- kDelta padrows) + // with polynom of 2nd order and two polynoms of 1st order + // take both polynoms of 1st order, calculate difference of their parameters + // add covariance matrixes and calculate chi2 of this difference + // if this chi2 is bigger than a given threshold, assume that the current cluster is + // a kink an goto next padrow + + + for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ + // loop over irow +- kDelta rows (neighboured rows) + // + // + if (idelta == 0) continue; + if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC + AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta); if (!cluster) continue; - if (cluster->GetType()<0) continue; - if (cluster->GetDetector()!=sector) continue; - Double_t x = cluster->GetX()-xref; + if (cluster->GetType() < 0) continue; + if (cluster->GetDetector() != sector) continue; + Double_t x = cluster->GetX() - xref; // x = differece: current cluster - cluster @ irow nclFound++; - if (idelta<0){ + if (idelta < 0){ ncl0++; - fitterY0.AddPoint(&x, cluster->GetY(),csigmaY); - fitterZ0.AddPoint(&x, cluster->GetZ(),csigmaZ); + fitterY0.AddPoint(&x, cluster->GetY(), csigmaY); + fitterZ0.AddPoint(&x, cluster->GetZ(), csigmaZ); } - if (idelta>0){ + if (idelta > 0){ ncl1++; - fitterY1.AddPoint(&x, cluster->GetY(),csigmaY); - fitterZ1.AddPoint(&x, cluster->GetZ(),csigmaZ); + fitterY1.AddPoint(&x, cluster->GetY(), csigmaY); + fitterZ1.AddPoint(&x, cluster->GetZ(), csigmaZ); } - fitterY2.AddPoint(&x, cluster->GetY(),csigmaY); - fitterZ2.AddPoint(&x, cluster->GetZ(),csigmaZ); - } - if (nclFoundGetY(), csigmaY); + fitterZ2.AddPoint(&x, cluster->GetZ(), csigmaZ); + } // loop over neighbourhood for fitter filling + + if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters found in neighbourhood goto next padrow fitterY2.Eval(); fitterZ2.Eval(); - Double_t chi2 = (fitterY2.GetChisquare()+fitterZ2.GetChisquare())/(2.*nclFound-6.); - if (chi2>kCutChi2) continue; - // - // + Double_t chi2 = (fitterY2.GetChisquare() + fitterZ2.GetChisquare()) / (2. * nclFound - 6.); + if (chi2 > kCutChi2) continue; // if chi^2 is too big goto next padrow + // REMOVE KINK - // - if (ncl0>4){ + // only when there are enough clusters (4) in each direction + if (ncl0 > 4){ fitterY0.Eval(); fitterZ0.Eval(); } - if (ncl1>4){ + if (ncl1 > 4){ fitterY1.Eval(); fitterZ1.Eval(); } - // - // - if (ncl0>4&&ncl1>4){ + + if (ncl0 > 4 && ncl1 > 4){ fitterY0.GetCovarianceMatrix(matrixY0); fitterY1.GetCovarianceMatrix(matrixY1); fitterZ0.GetCovarianceMatrix(matrixZ0); @@ -438,25 +640,37 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ fitterZ1.GetParameters(paramZ1); fitterY0.GetParameters(paramY0); fitterZ0.GetParameters(paramZ0); - paramY0-= paramY1; - paramZ0-= paramZ1; - matrixY0+= matrixY1; - matrixZ0+= matrixZ1; - Double_t chi2 =0; - TMatrixD difY(2,1,paramY0.GetMatrixArray()); - TMatrixD difYT(1,2,paramY0.GetMatrixArray()); + paramY0 -= paramY1; + paramZ0 -= paramZ1; + matrixY0 += matrixY1; + matrixZ0 += matrixZ1; + Double_t chi2 = 0; + + TMatrixD difY(2, 1, paramY0.GetMatrixArray()); + TMatrixD difYT(1, 2, paramY0.GetMatrixArray()); matrixY0.Invert(); - TMatrixD mulY(matrixY0,TMatrixD::kMult,difY); - TMatrixD chi2Y(difYT,TMatrixD::kMult,mulY); - chi2+=chi2Y(0,0); - TMatrixD difZ(2,1,paramZ0.GetMatrixArray()); - TMatrixD difZT(1,2,paramZ0.GetMatrixArray()); + TMatrixD mulY(matrixY0, TMatrixD::kMult, difY); + TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY); + chi2 += chi2Y(0, 0); + + TMatrixD difZ(2, 1, paramZ0.GetMatrixArray()); + TMatrixD difZT(1, 2, paramZ0.GetMatrixArray()); matrixZ0.Invert(); - TMatrixD mulZ(matrixZ0,TMatrixD::kMult,difZ); - TMatrixD chi2Z(difZT,TMatrixD::kMult,mulZ); - chi2+= chi2Z(0,0); - if (chi2*0.25>kCutChi2) continue; - } + TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ); + TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ); + chi2 += chi2Z(0, 0); + + // REMOVE KINK + if (chi2 * 0.25 > kCutChi2) continue; // if chi2 is too big goto next padrow + // fit tracklet with polynom of 2nd order and two polynoms of 1st order + // take both polynoms of 1st order, calculate difference of their parameters + // add covariance matrixes and calculate chi2 of this difference + // if this chi2 is bigger than a given threshold, assume that the current cluster is + // a kink an goto next padrow + } + // current padrow has no kink + + // get fit parameters from pol2 fit: Double_t paramY[4], paramZ[4]; paramY[0] = fitterY2.GetParameter(0); paramY[1] = fitterY2.GetParameter(1); @@ -464,77 +678,68 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ paramZ[0] = fitterZ2.GetParameter(0); paramZ[1] = fitterZ2.GetParameter(1); paramZ[2] = fitterZ2.GetParameter(2); - // - // - // + Double_t tracky = paramY[0]; Double_t trackz = paramZ[0]; - Float_t deltay = tracky-cluster0->GetY(); - Float_t deltaz = trackz-cluster0->GetZ(); - Float_t angley = paramY[1]-paramY[0]/xref; + Float_t deltay = tracky - cluster0->GetY(); + Float_t deltaz = trackz - cluster0->GetZ(); + Float_t angley = paramY[1] - paramY[0] / xref; Float_t anglez = paramZ[1]; - // - // + Float_t max = cluster0->GetMax(); TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector); if ( irow >= kFirstLargePad) max /= kLargePadSize; - profAmpRow->Fill(cluster0->GetRow(), max); + profAmpRow->Fill( (Double_t)cluster0->GetRow(), max ); TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector); hisAmp->Fill(max); - // - // - Int_t ipad=0; - if (cluster0->GetDetector()>=36) { - ipad=1; - if (cluster0->GetRow()>63) ipad=2; + + Int_t ipad = 0; + if (cluster0->GetDetector() >= 36) { + ipad = 1; + if (cluster0->GetRow() > 63) ipad = 2; } - // - // - TH3F * his3=0; + + TH3F * his3 = 0; his3 = (TH3F*)fRMSY->At(ipad); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley),TMath::Sqrt(cluster0->GetSigmaY2())); + if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) ); his3 = (TH3F*)fRMSZ->At(ipad); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez),TMath::Sqrt(cluster0->GetSigmaZ2())); - his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(),ipad)); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2())); - // - his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(),ipad)); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez),TMath::Sqrt(cluster0->GetSigmaZ2())); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) ); + + his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), ipad)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) ); + his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), ipad)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) ); - // // Fill resolution histograms - // - Bool_t useForResol= kTRUE; - if (fitterY2.GetParError(0)>kErrorFraction*csigmaY) useForResol=kFALSE; + Bool_t useForResol = kTRUE; + if (fitterY2.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE; if (useForResol){ fDeltaY->Fill(deltay); fDeltaZ->Fill(deltaz); his3 = (TH3F*)fResolY->At(ipad); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay ); his3 = (TH3F*)fResolZ->At(ipad); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz); - his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(),ipad)); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay); - // - his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(),ipad)); - if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz ); + his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), ipad)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay ); + his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), ipad)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz ); } - // - // - // - if (useForResol&&nclFound>2*kMinRatio*kDelta){ - // + + + if (useForResol && nclFound > 2 * kMinRatio * kDelta){ + // // fill resolution trees - // - static TLinearFitter fitY0(3,"pol2"); - static TLinearFitter fitZ0(3,"pol2"); - static TLinearFitter fitY2(5,"hyp4"); - static TLinearFitter fitZ2(5,"hyp4"); - static TLinearFitter fitY2Q(5,"hyp4"); - static TLinearFitter fitZ2Q(5,"hyp4"); - static TLinearFitter fitY2S(5,"hyp4"); - static TLinearFitter fitZ2S(5,"hyp4"); + // + static TLinearFitter fitY0(3, "pol2"); + static TLinearFitter fitZ0(3, "pol2"); + static TLinearFitter fitY2(5, "hyp4"); + static TLinearFitter fitZ2(5, "hyp4"); + static TLinearFitter fitY2Q(5, "hyp4"); + static TLinearFitter fitZ2Q(5, "hyp4"); + static TLinearFitter fitY2S(5, "hyp4"); + static TLinearFitter fitZ2S(5, "hyp4"); fitY0.ClearPoints(); fitZ0.ClearPoints(); fitY2.ClearPoints(); @@ -544,44 +749,45 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ fitY2S.ClearPoints(); fitZ2S.ClearPoints(); - for (Int_t idelta=-kDelta; idelta<=kDelta; idelta++){ - if (idelta==0) continue; - if (idelta+irow<0) continue; - if (idelta+irow>159) continue; - AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta); + for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ + // loop over irow +- kDelta rows (neighboured rows) + // + // + if (idelta == 0) continue; + if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC +// if (idelta + irow > 159) continue; + AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta); if (!cluster) continue; - if (cluster->GetType()<0) continue; - if (cluster->GetDetector()!=sector) continue; - Double_t x = cluster->GetX()-xref; - Double_t sigmaY0 = fClusterParam->GetError0Par(0,ipad,(250.0-TMath::Abs(cluster->GetZ())),TMath::Abs(angley)); - Double_t sigmaZ0 = fClusterParam->GetError0Par(1,ipad,(250.0-TMath::Abs(cluster->GetZ())),TMath::Abs(anglez)); + if (cluster->GetType() < 0) continue; + if (cluster->GetDetector() != sector) continue; + Double_t x = cluster->GetX() - xref; + Double_t sigmaY0 = fClusterParam->GetError0Par( 0, ipad, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) ); + Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, ipad, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) ); // - Double_t sigmaYQ = fClusterParam->GetErrorQPar(0,ipad,(250.0-TMath::Abs(cluster->GetZ())), - TMath::Abs(angley), TMath::Abs(cluster->GetMax())); - Double_t sigmaZQ = fClusterParam->GetErrorQPar(1,ipad,(250.0-TMath::Abs(cluster->GetZ())), - TMath::Abs(anglez),TMath::Abs(cluster->GetMax())); - Double_t sigmaYS = fClusterParam->GetErrorQParScaled(0,ipad,(250.0-TMath::Abs(cluster->GetZ())), - TMath::Abs(angley), TMath::Abs(cluster->GetMax())); - Double_t sigmaZS = fClusterParam->GetErrorQParScaled(1,ipad,(250.0-TMath::Abs(cluster->GetZ())), - TMath::Abs(anglez),TMath::Abs(cluster->GetMax())); - Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,ipad,(250.0-TMath::Abs(cluster->GetZ())), + Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, ipad, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) ); + Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, ipad, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) ); + Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, ipad, (250.0 - TMath::Abs(cluster->GetZ())), + TMath::Abs(angley), TMath::Abs(cluster->GetMax()) ); + Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, ipad, (250.0 - TMath::Abs(cluster->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) ); + Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, ipad,(250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()), - TMath::Sqrt(cluster0->GetSigmaY2()),0); - Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,ipad,(250.0-TMath::Abs(cluster->GetZ())), + TMath::Sqrt(cluster0->GetSigmaY2()), 0 ); + Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, ipad,(250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()), - TMath::Sqrt(cluster0->GetSigmaZ2()),0); - sigmaYS = TMath::Sqrt(sigmaYS*sigmaYS+rmsYFactor*rmsYFactor/12.); - sigmaZS = TMath::Sqrt(sigmaZS*sigmaZS+rmsZFactor*rmsZFactor/12.+rmsYFactor*rmsYFactor/24.); + TMath::Sqrt(cluster0->GetSigmaZ2()),0 ); + sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.); + sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.); // - if (kDelta!=0){ + if (kDelta != 0){ fitY0.AddPoint(&x, cluster->GetY(), sigmaY0); fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0); } Double_t xxx[4]; - xxx[0] = ((idelta+irow)%2==0)? 1:0; - xxx[1] = x; - xxx[2] = ((idelta+irow)%2==0)? x:0; - xxx[3] = x*x; + xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0; + xxx[1] = x; + xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0; + xxx[3] = x * x; fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0); fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ); fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS); @@ -589,7 +795,7 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ); fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS); // - } + } // neigbouhood-loop // fitY0.Eval(); fitZ0.Eval(); @@ -599,19 +805,19 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ fitZ2Q.Eval(); fitY2S.Eval(); fitZ2S.Eval(); - Float_t chi2Y0 = fitY0.GetChisquare()/(nclFound-3.); - Float_t chi2Z0 = fitZ0.GetChisquare()/(nclFound-3.); - Float_t chi2Y2 = fitY2.GetChisquare()/(nclFound-5.); - Float_t chi2Z2 = fitZ2.GetChisquare()/(nclFound-5.); - Float_t chi2Y2Q = fitY2Q.GetChisquare()/(nclFound-5.); - Float_t chi2Z2Q = fitZ2Q.GetChisquare()/(nclFound-5.); - Float_t chi2Y2S = fitY2S.GetChisquare()/(nclFound-5.); - Float_t chi2Z2S = fitZ2S.GetChisquare()/(nclFound-5.); + Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.); + Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.); + Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.); + Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.); + Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.); + Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.); + Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.); + Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.); // static TVectorD parY0(3); - static TMatrixD matY0(3,3); + static TMatrixD matY0(3, 3); static TVectorD parZ0(3); - static TMatrixD matZ0(3,3); + static TMatrixD matZ0(3, 3); fitY0.GetParameters(parY0); fitY0.GetCovarianceMatrix(matY0); fitZ0.GetParameters(parZ0); @@ -658,9 +864,8 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1)); Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3)); Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3)); - // + // Error parameters - // Float_t csigmaY0 = fClusterParam->GetError0Par(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley)); Float_t csigmaZ0 = fClusterParam->GetError0Par(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez)); // @@ -672,15 +877,15 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(anglez),TMath::Abs(cluster0->GetMax())); - // + // RMS parameters - // Float_t meanRMSY = 0; Float_t meanRMSZ = 0; - Int_t nclRMS=0; - for (Int_t idelta=-2; idelta<=2; idelta++){ - if (idelta+irow<0) continue; - if (idelta+irow>159) continue; + Int_t nclRMS = 0; + for (Int_t idelta = -2; idelta <= 2; idelta++){ + // loop over neighbourhood + if (idelta+irow < 0 || idelta+irow > 159) continue; +// if (idelta+irow>159) continue; AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta); if (!cluster) continue; meanRMSY += TMath::Sqrt(cluster->GetSigmaY2()); @@ -690,11 +895,11 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ meanRMSY /= nclRMS; meanRMSZ /= nclRMS; - Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2()); - Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2()); - Float_t rmsYT = fClusterParam->GetRMSQ(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())), + Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2()); + Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2()); + Float_t rmsYT = fClusterParam->GetRMSQ(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); - Float_t rmsZT = fClusterParam->GetRMSQ(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())), + Float_t rmsZT = fClusterParam->GetRMSQ(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster0->GetMax())); Float_t rmsYT0 = fClusterParam->GetRMS0(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(angley)); @@ -710,188 +915,187 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()), rmsZ,meanRMSZ); - // + // cluster debug - // - (*fDebugStream)<<"ResolCl"<< - "Sector="< 2 * kMinRatio * kDelta) + } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++) +} // FillResolutionHistoLocal(...) void AliTPCcalibTracks::AlignUpDown(AliTPCseed * track, AliESDtrack * esdTrack){ // // Make simple parabolic fit + // Write debug information to 'TPCSelectorDebug.root' // const Int_t kMinClusters = 60; - const Int_t kMinClustersSector =15; + const Int_t kMinClustersSector = 15; const Float_t kSigmaCut = 6; - const Float_t kMaxTan = TMath::Tan(TMath::Pi()*10./180.); + const Float_t kMaxTan = TMath::Tan(TMath::Pi() * 10. / 180.); const Float_t kDeadZone = 6.; const Float_t kMinZ = 15; - if (track->GetNumberOfClusters()GetZ())GetNumberOfClusters() < kMinClusters) return; + if (TMath::Abs(track->GetZ()) < kMinZ) return; // Int_t nclUp = 0; Int_t nclDown = 0; Int_t rSector =-1; - Float_t refX = (param.GetInnerRadiusUp()+param.GetOuterRadiusLow())*0.5; - for (Int_t irow=0; irow<159; irow++){ + Float_t refX = (param.GetInnerRadiusUp() + param.GetOuterRadiusLow()) * 0.5; + + for (Int_t irow = 0; irow < 159; irow++){ AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); if (!cluster0) continue; Int_t sector = cluster0->GetDetector(); - if (rSector<0) rSector=sector%36; - if (sector%36 != rSector) continue; - if ( ((TMath::Abs(cluster0->GetY())-kDeadZone)/cluster0->GetX())>kMaxTan) continue; //remove edge clusters - if (sector>35) nclUp++; - if (sector<36) nclDown++; - } - if (nclUpGetY()) - kDeadZone) / cluster0->GetX()) > kMaxTan ) continue; //remove edge clusters + if (sector > 35) nclUp++; + if (sector < 36) nclDown++; + } // loop over padrows + if (nclUp < kMinClustersSector) return; + if (nclDown < kMinClustersSector) return; + TLinearFitter fitterY(5,"hyp4"); //fitter with common 2 nd derivation TLinearFitter fitterZ(5,"hyp4"); - // - TLinearFitter fitterY0(3,"pol2"); + TLinearFitter fitterY0(3,"pol2"); // valgrind 3 58,142 bytes in 2,117 blocks are indirectly lost + // valgrind 608,151 (198,860 direct, 409,291 indirect) bytes in 1,108 blocks are definitely lost TLinearFitter fitterZ0(3,"pol2"); - TLinearFitter fitterY1(3,"pol2"); - TLinearFitter fitterZ1(3,"pol2"); - // - Float_t msigmay =1; - Float_t msigmaz =1; + TLinearFitter fitterY1(3,"pol2"); // valgrind 4,284 bytes in 21 blocks are possibly lost + TLinearFitter fitterZ1(3,"pol2"); // valgrind 6 blocks possibly lost 57,956 bytes in 844 blocks are indirectly lost + + Float_t msigmay = 1; + Float_t msigmaz = 1; Float_t param0[3]; Float_t param1[3]; - Float_t angley=0; - Float_t anglez=0; - // - for (Int_t iter=0; iter<3; iter++){ + Float_t angley = 0; + Float_t anglez = 0; + + for (Int_t iter = 0; iter < 3; iter++){ nclUp = 0; nclDown= 0; - for (Int_t irow=0; irow<159; irow++){ + for (Int_t irow = 0; irow < 159; irow++){ AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); if (!cluster0) continue; Int_t sector = cluster0->GetDetector(); - if (sector%36 != rSector) continue; + if (sector % 36 != rSector) continue; Double_t y = cluster0->GetY(); Double_t z = cluster0->GetZ(); //remove edge clusters - if ( (iter==0) && ((TMath::Abs(cluster0->GetY())-kDeadZone)/cluster0->GetX())>kMaxTan ) continue; - if (iter>0){ - Float_t tx = cluster0->GetX()-refX; + if ( (iter == 0) && ((TMath::Abs(cluster0->GetY()) - kDeadZone) / cluster0->GetX()) > kMaxTan ) continue; + if (iter > 0){ + Float_t tx = cluster0->GetX() - refX; Float_t ty = 0; - if (sector<36){ - ty = param0[0]+param0[1]*tx+param0[2]*tx*tx; + if (sector < 36){ + ty = param0[0] + param0[1] * tx + param0[2] * tx * tx; }else{ - ty = param1[0]+param1[1]*tx+param1[2]*tx*tx; + ty = param1[0] + param1[1] * tx + param1[2] * tx * tx; } if (((TMath::Abs(ty)-kDeadZone)/cluster0->GetX())>kMaxTan) continue; if (TMath::Abs(ty-y)>kSigmaCut*(msigmay+0.2)) continue; } - Int_t ipad=0; - if (cluster0->GetDetector()>=36) { - ipad=1; + Int_t ipad = 0; + if (cluster0->GetDetector() >= 36) { + ipad = 1; if (cluster0->GetRow()>63) ipad=2; } // - Float_t sigmaY =msigmay; - Float_t sigmaZ =msigmay; - if (iter==2){ + Float_t sigmaY = msigmay; + Float_t sigmaZ = msigmay; + if (iter == 2){ sigmaY = fClusterParam->GetErrorQParScaled(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); sigmaZ = fClusterParam->GetErrorQParScaled(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(anglez),TMath::Abs(cluster0->GetMax())); - } - Double_t deltaX = cluster0->GetX()-refX; + } // iter == 2 + Double_t deltaX = cluster0->GetX() - refX; Double_t x[5]; x[0] = (ipad==0) ? 0:1; x[1] = deltaX; x[2] = (ipad==0) ? 0:deltaX; x[3] = deltaX*deltaX; - if (ipad<2){ + if (ipad < 2){ fitterY.AddPoint(x,y,sigmaY); fitterZ.AddPoint(x,z,sigmaZ); } - if (ipad==0){ + if (ipad == 0){ nclDown++; fitterY0.AddPoint(&deltaX,y,sigmaY); fitterZ0.AddPoint(&deltaX,z,sigmaZ); } - if (ipad==1){ + if (ipad == 1){ nclUp++; fitterY1.AddPoint(&deltaX,y,sigmaY); fitterZ1.AddPoint(&deltaX,z,sigmaZ); } - } - if (nclUpGetFriendTrack(); AliTrackPointArray *points = (AliTrackPointArray*)ftrack->GetTrackPointArray(); - - if (iter>0) (*fDebugStream)<<"Align"<< + + if (iter>0) (*fDebugStream)<<"Align"<< // valgrind 85,932 bytes in 543 blocks are still reachable "track.="<Clone(); // valgrind 3 40,139 bytes in 11 blocks are still reachable + result->SetTitle(Form("%s fit residuals",result->GetTitle())); + result->SetName(Form("%s fit residuals",result->GetName())); + TAxis *xaxis = hfit->GetXaxis(); + TAxis *yaxis = hfit->GetYaxis(); + Double_t x[2]; + for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) { + x[1] = yaxis->GetBinCenter(biny); + for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) { + x[0] = xaxis->GetBinCenter(binx); + Int_t bin = hfit->GetBin(binx, biny); + Double_t val = hfit->GetBinContent(bin); +// result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) ); + result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 ); + } + } + return result; +} + + + +void AliTPCcalibTracks::SetStyle(){ + // + // set style, can be called by all draw functions + // + gROOT->SetStyle("Plain"); + gStyle->SetFillColor(10); + gStyle->SetPadColor(10); + gStyle->SetCanvasColor(10); + gStyle->SetStatColor(10); + gStyle->SetPalette(1,0); + gStyle->SetNumberContours(60); +} + + +void AliTPCcalibTracks::Draw(Option_t* opt){ + // + // draw-function of AliTPCcalibTracks + // will draws some exemplaric pictures + // + + cout << "***** not yet implemented *****" << endl; + fDeltaY->Draw(opt); + +} + + +void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ + // + // all functions are called, that produce pictures + // the histograms are written to the directory 'pathName' + // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file + // 'stat' is also the number of minEntries for MakeResPlotsQTree + // + + + MakeAmpPlots(stat, pathName); + MakeDeltaPlots(pathName); + FitResolutionNew(pathName); + FitRMSNew(pathName); +// cout << "now comes MakeResPlotsQ" << endl; +// MakeResPlotsQ(1, 1); +// cout << "now comes MakeResPlotsQTree" << endl; + MakeResPlotsQTree(stat, pathName); + +} + + +void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ + // + // creates several plots: + // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps + // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster + // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit + // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit + // Empty histograms (sectors without data) are not written to file + // the ps-files are written to the directory 'pathName', that is created if it does not exist + // 'stat': only histograms with more than 'stat' entries are written to file. + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable + TPostScript *ps; + // histograms with accumulated amplitude for all IROCs and OROCs + TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone()); + allAmpHisIROC->SetName("Amp all IROCs"); + allAmpHisIROC->SetTitle("Amp all IROCs"); + TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone()); + allAmpHisOROC->SetName("Amp all OROCs"); + allAmpHisOROC->SetTitle("Amp all OROCs"); + + + ps = new TPostScript("fArrayAmp.ps", 112); + cout << "creating fArrayAmp.ps..." << endl; + for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){ + if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat ) continue; + ps->NewPage(); + ((TH1F*)fArrayAmp->At(i))->Draw(); + c1->Update(); // valgrind 3 + if (i > 0 && i < 36) { + allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i))); + allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36))); + } + } + ps->NewPage(); + allAmpHisIROC->Draw(); + c1->Update(); // valgrind + ps->NewPage(); + allAmpHisOROC->Draw(); + c1->Update(); + ps->Close(); + delete ps; + + TH1F *his = 0; + Double_t min = 0; + Double_t max = 0; + ps = new TPostScript("fArrayAmpRow.ps", 112); + cout << "creating fArrayAmpRow.ps..." << endl; + for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){ + his = (TH1F*)fArrayAmpRow->At(i); + if (his->GetEntries() < stat) continue; + ps->NewPage(); + min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.); + max = his->GetBinCenter(5*his->GetMaximumBin()) + 100; + his->SetAxisRange(min, max); + his->Fit("pol3", "q", "", min, max); + // his->Draw("error"); // don't use this line when you don't want to have empty pages in the ps-file + c1->Update(); + } + ps->Close(); + delete ps; + + delete c1; + gSystem->ChangeDirectory(".."); +} + + +void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){ + // + // creates several plots: + // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit + // the ps-files are written to the directory 'pathName', that is created if it does not exist + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable + TPostScript *ps; + Double_t min = 0; + Double_t max = 0; + + ps = new TPostScript("DeltaYZ.ps", 112); + cout << "creating DeltaYZ.ps..." << endl; + min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20; + max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20; + fDeltaY->SetAxisRange(min, max); + ps->NewPage(); +// fDeltaY->Draw(); + fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable + c1->Update(); + ps->NewPage(); + max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20; + min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20; + fDeltaZ->SetAxisRange(min, max); +// fDeltaZ->Draw(); + fDeltaZ->Fit("gaus","q","",min, max); + c1->Update(); + ps->Close(); + delete ps; + delete c1; + + gSystem->ChangeDirectory(".."); +} + + +void AliTPCcalibTracks::FitResolutionNew(char* pathName){ + // + // calculates different resulution fits in Y and Z direction + // the histograms are written to 'ResolutionYZ.ps' + // writes calculated resolution to 'resol.txt' + // all files are stored in the directory pathName + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas c; + c.Divide(2,1); + cout << "creating ResolutionYZ.ps..." << endl; + TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112); + TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1); + fres->SetParameter(0,0.02); + fres->SetParameter(1,0.0054); + fres->SetParameter(2,0.13); + + TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory + + // create histogramw for Y-resolution + TH3F * hisResY0 = (TH3F*)fResolY->At(0); + hisResY0->FitSlicesZ(); + TH2D * hisResY0_2 = (TH2D*)gDirectory->Get("Resol Y0_2"); + TH3F * hisResY1 = (TH3F*)fResolY->At(1); + hisResY1->FitSlicesZ(); + TH2D * hisResY1_2 = (TH2D*)gDirectory->Get("Resol Y1_2"); + TH3F * hisResY2 = (TH3F*)fResolY->At(2); + hisResY2->FitSlicesZ(); + TH2D * hisResY2_2 = (TH2D*)gDirectory->Get("Resol Y2_2"); + // + ps->NewPage(); + c.cd(1); + hisResY0_2->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost + hisResY0_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY0_2,fres)->Draw("surf1"); + c.Update(); + // c.SaveAs("ResolutionYPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResY1_2->Fit(fres, "q"); + hisResY1_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY1_2,fres)->Draw("surf1"); + c.Update(); + // c.SaveAs("ResolutionYPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResY2_2->Fit(fres, "q"); + hisResY2_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY2_2,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionYPad2.eps"); + + // create histogramw for Z-resolution + TH3F * hisResZ0 = (TH3F*)fResolZ->At(0); + hisResZ0->FitSlicesZ(); + TH2D * hisResZ0_2 = (TH2D*)gDirectory->Get("Resol Z0_2"); + TH3F * hisResZ1 = (TH3F*)fResolZ->At(1); + hisResZ1->FitSlicesZ(); + TH2D * hisResZ1_2 = (TH2D*)gDirectory->Get("Resol Z1_2"); + TH3F * hisResZ2 = (TH3F*)fResolZ->At(2); + hisResZ2->FitSlicesZ(); + TH2D * hisResZ2_2 = (TH2D*)gDirectory->Get("Resol Z2_2"); + + ps->NewPage(); + c.cd(1); + hisResZ0_2->Fit(fres, "q"); + hisResZ0_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ0_2,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionZPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResZ1_2->Fit(fres, "q"); + hisResZ1_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ1_2,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionZPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResZ2_2->Fit(fres, "q"); + hisResZ2_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ2_2,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionZPad2.eps"); + ps->Close(); + delete ps; + + // write calculated resoltuions to 'resol.txt' + ofstream fresol("resol.txt"); + fresol<<"Pad 0.75 cm"<<"\n"; + hisResY0_2->Fit(fres, "q"); // valgrind + fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ0_2->Fit(fres, "q"); + fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + fresol<<"Pad 1.00 cm"<<1<<"\n"; + hisResY1_2->Fit(fres, "q"); // valgrind + fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ1_2->Fit(fres, "q"); + fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + fresol<<"Pad 1.50 cm"<<0<<"\n"; + hisResY2_2->Fit(fres, "q"); + fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ2_2->Fit(fres, "q"); + fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + + TH1::AddDirectory(kFALSE); + gSystem->ChangeDirectory(".."); + delete fres; +} + + +void AliTPCcalibTracks::FitRMSNew(char* pathName){ + // + // calculates different resulution-rms fits in Y and Z direction + // the histograms are written to 'RMS_YZ.ps' + // writes calculated resolution rms to 'rms.txt' + // all files are stored in the directory pathName + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas c; // valgrind 3 42,120 bytes in 405 blocks are still reachable 23,816 bytes in 229 blocks are still reachable + c.Divide(2,1); + cout << "creating RMS_YZ.ps..." << endl; + TPostScript *ps = new TPostScript("RMS_YZ.ps", 112); + TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1); + frms->SetParameter(0,0.02); + frms->SetParameter(1,0.0054); + frms->SetParameter(2,0.13); + + TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory + + // create histogramw for Y-RMS + TH3F * hisResY0 = (TH3F*)fRMSY->At(0); + hisResY0->FitSlicesZ(); + TH2D * hisResY0_2 = (TH2D*)gDirectory->Get("RMS Y0_1"); + TH3F * hisResY1 = (TH3F*)fRMSY->At(1); + hisResY1->FitSlicesZ(); + TH2D * hisResY1_2 = (TH2D*)gDirectory->Get("RMS Y1_1"); + TH3F * hisResY2 = (TH3F*)fRMSY->At(2); + hisResY2->FitSlicesZ(); + TH2D * hisResY2_2 = (TH2D*)gDirectory->Get("RMS Y2_1"); + // + ps->NewPage(); + c.cd(1); + hisResY0_2->Fit(frms, "qn0"); + hisResY0_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY0_2,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSYPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResY1_2->Fit(frms, "qn0"); // valgrind several blocks possibly lost + hisResY1_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY1_2,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSYPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResY2_2->Fit(frms, "qn0"); + hisResY2_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY2_2,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSYPad2.eps"); + + // create histogramw for Z-RMS + TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0); + hisResZ0->FitSlicesZ(); + TH2D * hisResZ0_2 = (TH2D*)gDirectory->Get("RMS Z0_1"); + TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1); + hisResZ1->FitSlicesZ(); + TH2D * hisResZ1_2 = (TH2D*)gDirectory->Get("RMS Z1_1"); + TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2); + hisResZ2->FitSlicesZ(); + TH2D * hisResZ2_2 = (TH2D*)gDirectory->Get("RMS Z2_1"); + // + ps->NewPage(); + c.cd(1); + hisResZ0_2->Fit(frms, "qn0"); // valgrind + hisResZ0_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ0_2,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSZPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResZ1_2->Fit(frms, "qn0"); + hisResZ1_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ1_2,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSZPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResZ2_2->Fit(frms, "qn0"); // valgrind 1 block possibly lost + hisResZ2_2->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ2_2,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSZPad2.eps"); + + // write calculated resoltuion rms to 'rms.txt' + ofstream filerms("rms.txt"); + filerms<<"Pad 0.75 cm"<<"\n"; + hisResY0_2->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost + filerms<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ0_2->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost + filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + filerms<<"Pad 1.00 cm"<<1<<"\n"; + hisResY1_2->Fit(frms, "qn0"); // valgrind 3,256 bytes in 22 blocks are indirectly lost + filerms<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ1_2->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable + filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + filerms<<"Pad 1.50 cm"<<0<<"\n"; + hisResY2_2->Fit(frms, "qn0"); // valgrind 40,139 bytes in 11 blocks are still reachable 330,180 bytes in 15 blocks are possibly lost + filerms<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ2_2->Fit(frms, "qn0"); + filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + + TH1::AddDirectory(kFALSE); + gSystem->ChangeDirectory(".."); + ps->Close(); + delete ps; + delete frms; +} + + + +void AliTPCcalibTracks::MakeResPlotsQ(Int_t minEntries, Bool_t bDraw, char* pathName){ + // + // make resolution Plots + // not yet finished function + // + cout << " not yet finished function MakeResPlotsQ" << endl; + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TLinearFitter fitter(3,"hyp2"); + TLinearFitter fitterQ(3,"hyp2"); + Int_t npointsFit =0; + TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1); + fres->SetParameter(0,0.02); + fres->SetParameter(1,0.0054); + fres->SetParameter(2,0.13); + fres->SetParLimits(0,0,0.1); + TGraph2DErrors *graph = new TGraph2DErrors(1000); + + for (Int_t idim = 0; idim < 2; idim++){ + const char* dim = (idim==0)? "Y":"Z"; + for (Int_t ipad = 0; ipad < 3; ipad++){ + // + printf("Direction %s\t",dim); + printf("Pad %d\n",ipad); + npointsFit = 0; + fitter.ClearPoints(); + fitterQ.ClearPoints(); + for (Int_t iq = 0; iq < 4; iq++){ + Int_t bin = GetBin(iq, ipad); + Float_t qmean = GetQ(bin); + char name[200]; + sprintf(name, "Resol%s Pad%d Qmiddle%f",dim, ipad, qmean); + TH3F *hisin = 0; + if (idim == 0) hisin = (TH3F*)fArrayQDY->At(bin); + if (idim == 1) hisin = (TH3F*)fArrayQDZ->At(bin); + if (!hisin) continue; + //printf("Q\t%f\t",qmean); + cout << "calling FitProjections" << endl; + TObjArray * array = FitProjections(hisin, 2, minEntries, bDraw); + // + // Fit resolution + // + cout << "fit resolution" << endl; + cout << "array->GetEntriesFast(): " << array->GetEntriesFast() << endl; + //array->Print(); + for (Int_t ipoint = 0; ipoint < array->GetEntriesFast(); ipoint++){ + TVectorF * vector = (TVectorF*)array->At(ipoint); + if (!vector) continue; + //vector->Print(); + Double_t val = (*vector)[4]; + Double_t error = (*vector)[5]; + error *= 2 * val; + //error += val * val * 0.05; // 5% error addition + val *= val; + // + // mean fitter + // + Double_t zmean = (*vector)[0]; + Double_t angle = (*vector)[2]* (*vector)[2]; + Double_t x[2]; + x[0] = zmean; + x[1] = angle; + fitter.AddPoint(x,val,error); + // + // with q fitter + // + Double_t zmeanq = (*vector)[0]/qmean; + Double_t angleq = (*vector)[2]* (*vector)[2]; + Double_t xq[2]; + xq[0] = zmeanq; + xq[1] = angleq; + fitterQ.AddPoint(xq,val,error); + graph->SetPoint(npointsFit,(*vector)[0],(*vector)[2],(*vector)[4]); + graph->SetPointError(npointsFit,(*vector)[1],(*vector)[3],(*vector)[5]); + npointsFit++; + } + } + printf("NPoints = %d \n", npointsFit); + cout << "evaluating fitters" << endl; + fitter.Eval(); + fitterQ.Eval(); + TGraph2DErrors *graphVal = new TGraph2DErrors(npointsFit); + TGraph2DErrors *graphDif = new TGraph2DErrors(npointsFit); +/* graphVal->SetDirectory(0); + graphDif->SetDirectory(0);*/ + for (Int_t ipoint=0; ipointGetX()[ipoint]+0.01*Float_t(ipad); + x[1] = graph->GetY()[ipoint]+0.01*ipad; + x[2] = graph->GetZ()[ipoint]; + sigma[0] = graph->GetErrorX(ipoint); + sigma[1] = graph->GetErrorY(ipoint); + sigma[3] = graph->GetErrorZ(ipoint); + graphVal->SetPoint(ipoint,x[0],x[1],x[2]); + graphVal->SetPointError(ipoint,sigma[0],sigma[1],sigma[2]); + } + fres->SetParameter(0,0.02); + fres->SetParameter(1,0.0054); + fres->SetParameter(2,0.13); + + graphVal->Fit(fres,"q"); + for (Int_t ipoint=0; ipointGetX()[ipoint]+0.01*Float_t(ipad); + x[1] = graph->GetY()[ipoint]+0.01*ipad; + x[2] = graph->GetZ()[ipoint]-fres->Eval(x[0],x[1]); + sigma[0] = graph->GetErrorX(ipoint); + sigma[1] = graph->GetErrorY(ipoint); + sigma[3] = graph->GetErrorZ(ipoint); + graphDif->SetPoint(ipoint,x[0],x[1],x[2]); + graphDif->SetPointError(ipoint,sigma[0],sigma[1],sigma[2]); + } + Double_t p0 = TMath::Sqrt(TMath::Abs(fitter.GetParameter(0))); + Double_t p1 = TMath::Sqrt(TMath::Abs(fitter.GetParameter(1))); + Double_t p2 = TMath::Sqrt(TMath::Abs(fitter.GetParameter(2))); + Double_t chi2= fitter.GetChisquare()/npointsFit; + printf("Linear fit - chi2 %f\t%f\t%f\t%f\n",chi2,p0,p1,p2); + Double_t p0q = TMath::Sqrt(TMath::Abs(fitterQ.GetParameter(0))); + Double_t p1q = TMath::Sqrt(TMath::Abs(fitterQ.GetParameter(1))); + Double_t p2q = TMath::Sqrt(TMath::Abs(fitterQ.GetParameter(2))); + Double_t chi2q= fitterQ.GetChisquare()/npointsFit; + printf("Linear fitQ - chi2 %f\t%f\t%f\t%f\n",chi2q,p0q,p1q,p2q); + + printf("Graph fit - chi2 %f\t%f\t%f\t%f\n",fres->GetChisquare(),fres->GetParameter(0),fres->GetParameter(1),fres->GetParameter(2)); + } + } + gSystem->ChangeDirectory(".."); +} + + +TObjArray* AliTPCcalibTracks::FitProjections(TH3F * hfit, Int_t val, Int_t minEntry, Bool_t bDraw){ + // + // ??? + // needed by MakeResPlotsQ + // + + cout << "AliTPCcalibTracks::FitProjections started" << endl; + TObjArray * array = new TObjArray(100); + TAxis *xaxis = hfit->GetXaxis(); + TAxis *yaxis = hfit->GetYaxis(); + Double_t x[2]; + char name[200]; + sprintf(name,"Histos%s_%d", hfit->GetName(),val); + TCanvas *canvas = 0; + if (bDraw){ + canvas = new TCanvas(name,name); + canvas->Divide(yaxis->GetNbins(),xaxis->GetNbins()); } + Int_t count =1; + Int_t zaehler = 0; + for (Int_t biny = 1; biny <= yaxis->GetNbins(); biny++) { + x[1] = yaxis->GetBinCenter(biny); + for (Int_t binx = 1; binx <= xaxis->GetNbins(); binx++) { + x[0] = xaxis->GetBinCenter(binx); + char name[200]; + sprintf(name, "%s x %f y %f", hfit->GetName(),x[0],x[1]); + TH1D * projection = (TH1D*)( hfit->ProjectionZ(name, binx, binx, biny, biny) ); + projection->SetDirectory(0); + // + // + if (projection->GetEntries() < minEntry){ + Float_t meanz = 0; + Float_t meanangle = 0; + Double_t entries = 0; + projection->Clear(); + for (Int_t dbin = 0; dbin <= 4; dbin++) + for (Int_t dbiny2 = -3; dbiny2 <= 3; dbiny2++) { + for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ + zaehler++; + if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; + Int_t binx2 = binx + dbinx2; + Int_t biny2 = biny + dbiny2; +// if (binx2 < 1 || biny2 < 1 || binx2 > xaxis->GetNbins() || biny2 > yaxis->GetNbins()) continue; + if (binx2 < 1) continue; + if (biny2 < 1) continue; + if (binx2 > xaxis->GetNbins()) continue; + if (biny2 > yaxis->GetNbins()) continue; + TH1D * projection2 = (TH1D*)(hfit->ProjectionZ("Temp", binx2, binx2, biny2, biny2)); + //projection2->SetDirectory(0); + projection->Add(projection2); + entries += projection2->GetEntries(); + meanz += projection2->GetEntries() * xaxis->GetBinCenter(binx2); + meanangle += projection2->GetEntries() * yaxis->GetBinCenter(biny2); + if (entries > minEntry) break; + } + if (entries > minEntry) break; + } + if ( projection->GetEntries()GetEntries() < minEntry) continue; + // + Float_t xmin = projection->GetMean()-2.*projection->GetRMS()-0.02; + Float_t xmax = projection->GetMean()+2.*projection->GetRMS()+0.02; + // printf("%f\t%f\n",xmin,xmax); + if (bDraw) canvas->cd(count); + projection->Fit("gaus","q","",xmin,xmax); + TVectorF *pvector = new TVectorF(6); + TVectorF & vector = *pvector; + vector[0] = x[0]; + vector[1] = xaxis->GetBinWidth(binx); + vector[2] = x[1]; + vector[3] = yaxis->GetBinWidth(biny); + vector[4] = projection->GetFunction("gaus")->GetParameter(2); + vector[5] = projection->GetFunction("gaus")->GetParError(2); + array->AddLast(pvector); + count++; + } + } + +// cout << "the inner loop was executed "<< zaehler << " times!" << endl; + // + // + // + + TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory + +// cout << "critical point in FitProjections" << endl; + + TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1); + fres->SetParameter(0,0.02); + fres->SetParameter(1,0.0054); + fres->SetParameter(2,0.13); + // + hfit->FitSlicesZ(); + sprintf(name,"%s_%d", hfit->GetName(),val); + TH2D * his_2 = (TH2D*)gDirectory->Get(name); + his_2->Fit(fres,"q"); + his_2->SetDirectory(0); + if (bDraw){ + TCanvas *canvas2 = new TCanvas( hfit->GetName(), hfit->GetName()); + canvas2->Divide(1,2); + canvas2->cd(1); + his_2->Draw("surf1"); + canvas2->cd(2); + MakeDiff(his_2,fres)->Draw("surf1"); + } + // + printf("fres_Parameter[0]: %f\tfres_Parameter[1]: %f\tfres_Parameter[2]: %f\n",fres->GetParameter(0),fres->GetParameter(1),fres->GetParameter(2)); + + TH1::AddDirectory(kFALSE); + cout << "end of FitProjections, array->GetEntriesFast():" << array->GetEntriesFast() << endl; + return array; +} + + + +void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ + // + // Make tree with resolution parameters + // the result is written to 'resol.root' in directory 'pathname' + // + + cout << " ***** this is MakeResPlotsQTree *****" << endl; + cout << " relax, the calculation will take a while..." << endl; + + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + TTreeSRedirector fTreeResol("resol.root"); + + TH3F *resArray[2][3][11]; + TH3F *rmsArray[2][3][11]; + + // load histograms into resArraz and rmsArray + for (Int_t idim = 0; idim < 2; idim++){ + for (Int_t ipad = 0; ipad < 3; ipad++){ + for (Int_t iq = 0; iq <= 10; iq++){ + rmsArray[idim][ipad][iq]=0; + resArray[idim][ipad][iq]=0; + Int_t bin = GetBin(iq,ipad); +// Double_t qCenter = GetQ(bin); // unused variable !!! + + TH3F *hresl = 0; + if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin); + if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin); + if (!hresl) continue; + resArray[idim][ipad][iq] = (TH3F*) hresl->Clone(); + resArray[idim][ipad][iq]->SetDirectory(0); + + TH3F * hreslRMS = 0; + if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin); + if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin); + if (!hreslRMS) continue; + rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone(); + rmsArray[idim][ipad][iq]->SetDirectory(0); + } + } + } + + cout << "Histograms loaded, starting to proces..." << endl; + + //-------------------------------------------------------------------------------------------- + + //Double_t qCenter; + Double_t qMean; + Double_t zMean, angleMean, zCenter, angleCenter; + Double_t zSigma, angleSigma; + Int_t loopCounter = 1; + + for (Int_t idim = 0; idim < 2; idim++){ + // Loop y-z corrdinate + for (Int_t ipad = 0; ipad < 3; ipad++){ + // loop pad type + + // printf("%d\t%d\n", idim, ipad); + for (Int_t iq = -1; iq < 10; iq++){ + // LOOP Q + cout << "Loop-counter, this is loop " << loopCounter << " of 66, (" + << (Int_t)((loopCounter)/66.*100) << "% done), " + << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush; + loopCounter++; + + TH3F *hres = 0; + TH3F *hrms = 0; + qMean = 0; + if (iq == -1){ + // integrated spectra + Float_t entriesQ = 0; + for (Int_t iql = 0; iql < 10; iql++){ + Int_t bin = GetBin(iql,ipad); + TH3F *hresl = resArray[idim][ipad][iql]; + TH3F *hrmsl = rmsArray[idim][ipad][iql]; + if (!hresl) continue; + if (!hrmsl) continue; + entriesQ += hresl->GetEntries(); + qMean += hresl->GetEntries() * GetQ(bin); + if (!hres) { + hres = (TH3F*)hresl->Clone(); + hrms = (TH3F*)hrmsl->Clone(); + } + else{ + hres->Add(hresl); + hrms->Add(hrmsl); + } + } + qMean /= entriesQ; + qMean *= -1.; // integral mean charge + } + else { + Float_t entriesQ = 0; + for (Int_t iql = iq - 1; iql <= iq + 1; iql++){ + if (iql < 0) continue; + Int_t bin = GetBin(iql,ipad); + // qCenter = GetQ(bin); + TH3F * hresl = resArray[idim][ipad][iql]; + TH3F * hrmsl = rmsArray[idim][ipad][iql]; + if (!hresl) continue; + if (!hrmsl) continue; + entriesQ += hresl->GetEntries(); + qMean += hresl->GetEntries() * GetQ(bin); + if (!hres) { + hres = (TH3F*) hresl->Clone(); + hrms = (TH3F*) hrmsl->Clone(); + } + else{ + hres->Add(hresl); + hrms->Add(hrmsl); + } + } + qMean/=entriesQ; + } + + TAxis *xaxis = hres->GetXaxis(); + TAxis *yaxis = hres->GetYaxis(); + TAxis *zaxis = hres->GetZaxis(); + for (Int_t biny = 1; biny <= yaxis->GetNbins(); biny++) { + // angle loop + angleCenter = yaxis->GetBinCenter(biny); + for (Int_t binx = 1; binx <= xaxis->GetNbins(); binx++) { + // z - loop + zCenter = xaxis->GetBinCenter(binx); + zMean = zCenter; + angleMean = angleCenter; + zSigma = xaxis->GetBinWidth(binx); + angleSigma = yaxis->GetBinWidth(biny); + + + // create 2 1D-Histograms, projectionRes and projectionRms + // here it is possible to speed up the program by using the loop and the other + // TH1D *projectionRes... and TH1D *projectionRms... statements + // but there is a bug somewhere.... + char name[200]; + sprintf(name,"%s x %f y %f", hres->GetName(),zCenter,angleCenter); +// TH1D *projectionRes = new TH1D(name, name, zaxis->GetNbins(), zaxis->GetXmin(), zaxis->GetXmax()); + TH1D * projectionRes = (TH1D*)(hres->ProjectionZ(name,binx,binx, biny,biny)); + + sprintf(name,"%s x %f y %f", hrms->GetName(),zCenter,angleCenter); +// TH1D *projectionRms = new TH1D(name, name, zaxis->GetNbins(), zaxis->GetXmin(), zaxis->GetXmax()); + TH1D * projectionRms = (TH1D*)(hrms->ProjectionZ(name,binx,binx, biny,biny)); + +/* + for (Int_t ibin3 = 1; ibin3 < zaxis->GetNbins(); ibin3++) { +// fill 1D-Histograms + projectionRes->Fill(ibin3, hres->GetBinContent(binx, biny, ibin3)); + projectionRms->Fill(ibin3, hrms->GetBinContent(binx, biny, ibin3)); + } +*/ + projectionRes->SetDirectory(0); + projectionRms->SetDirectory(0); + + if (projectionRes->GetEntries() < minEntries){ + // if not enough statistic + zMean = 0; + angleMean = 0; + Double_t entries =0; + projectionRes->Clear(); + projectionRms->Clear(); + for (Int_t dbin = 0; dbin <= 8; dbin++) + for (Int_t dbiny2 = -4; dbiny2 <= 4; dbiny2++) { + for (Int_t dbinx2 = -4; dbinx2 <= 4; dbinx2++){ + if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; + Int_t binx2 = binx + dbinx2; + Int_t biny2 = biny + dbiny2; + if (binx2 < 1) continue; + if (biny2 < 1) continue; + if (binx2 >= xaxis->GetNbins()) continue; + if (biny2 >= yaxis->GetNbins()) continue; + // Code below is to slow +/* TH1D * projectionRes2 = (TH1D*)(hres->ProjectionZ("Temp",binx2,binx2, biny2,biny2)); + projectionRes2->SetDirectory(0); + projectionRes->Add(projectionRes2); + TH1D * projectionRMS2 = (TH1D*)(hrms->ProjectionZ("Temp",binx2,binx2, biny2,biny2)); + projectionRMS2->SetDirectory(0); + projectionRms->Add(projectionRMS2); + // + entries += projectionRes2->GetEntries(); + zMean += projectionRes2->GetEntries() * xaxis->GetBinCenter(binx2); + angleMean += projectionRes2->GetEntries() * yaxis->GetBinCenter(biny2); + delete projectionRes2; + delete projectionRMS2;*/ + for (Int_t ibin3 = 1; ibin3 < zaxis->GetNbins(); ibin3++) { + projectionRes->Fill(ibin3, hres->GetBinContent(binx2, biny2, ibin3)); + projectionRms->Fill(ibin3, hres->GetBinContent(binx2, biny2, ibin3)); + entries += hres->GetBinContent(binx2, biny2, ibin3); + zMean += hres->GetBinContent(binx2, biny2, ibin3) * xaxis->GetBinCenter(binx2); + angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yaxis->GetBinCenter(biny2); + } // ibin3 loop + if (entries > minEntries) break; + } //dbinx2 loop + if (entries > minEntries) break; + } // dbiny2 loop + if ( projectionRes->GetEntries() < minEntries) continue; + if (entries == 0) continue; + zMean /= entries; + angleMean /= entries; + } // if (projectionRes->GetEntries() < minEntries) + + if (projectionRes->GetEntries() > minEntries) { + // when enough statistic is accumulated + Float_t entries = projectionRes->GetEntries(); + Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.02; + Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.02; + projectionRes->Fit("gaus","q","",xmin,xmax); + Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2); + Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2); + // + xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.02; + xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.02; + projectionRms->Fit("gaus","q","",xmin,xmax); + Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1); + Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1); + Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2); + Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2); + // + Float_t length = 0.75; + if (ipad == 1) length = 1; + if (ipad == 2) length = 1.5; + + fTreeResol<<"Resol"<< + "Entries="<ChangeDirectory(".."); + } + + + +Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { + // + // function to merge several AliTPCcalibTracks objects after PROOF calculation + // The object's histograms are merged via their merge functions + // + + cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl; + if (!collectionList) return 0; + if (collectionList->IsEmpty()) return -1; + + cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1 + + // create a list for each data member + TList* deltaYList = new TList; + TList* deltaZList = new TList; + TList* arrayAmpRowList = new TList; + TList* arrayAmpList = new TList; + TList* arrayQDYList = new TList; + TList* arrayQDZList = new TList; + TList* arrayQRMSYList = new TList; + TList* arrayQRMSZList = new TList; + TList* resolYList = new TList; + TList* resolZList = new TList; + TList* rMSYList = new TList; + TList* rMSZList = new TList; + + TList* nRowsList = new TList; + TList* nSectList = new TList; + TList* fileNoList = new TList; + + TIterator *listIterator = collectionList->MakeIterator(); + AliTPCcalibTracks *calibTracks = 0; + cout << "start to iterate, filling lists" << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1 + Int_t counter = 0; + while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){ + // loop over all entries in the collectionList and get dataMembers into lists + if (!calibTracks) continue; + deltaYList->Add( calibTracks->GetfDeltaY() ); + deltaZList->Add( calibTracks->GetfDeltaZ() ); + arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow()); + arrayAmpList->Add(calibTracks->GetfArrayAmp()); + arrayQDYList->Add(calibTracks->GetfArrayQDY()); + arrayQDZList->Add(calibTracks->GetfArrayQDZ()); + arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY()); + arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ()); + resolYList->Add(calibTracks->GetfResolY()); + resolZList->Add(calibTracks->GetfResolZ()); + rMSYList->Add(calibTracks->GetfRMSY()); + rMSZList->Add(calibTracks->GetfRMSZ()); + counter++; + } + + // reset data members + cout << "histogram's reset-functins are called... " << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1 + fDeltaY->Reset(); + fDeltaZ->Reset(); + for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) + ((TProfile*)(fArrayAmpRow->At(i)))->Reset(); + for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) + ((TProfile*)(fArrayAmp->At(i)))->Reset(); + for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) + ((TH3F*)(fArrayQDY->At(i)))->Reset(); + for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) + ((TH3F*)(fArrayQDZ->At(i)))->Reset(); + for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) + ((TH3F*)(fArrayQRMSY->At(i)))->Reset(); + for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) + ((TH3F*)(fArrayQRMSZ->At(i)))->Reset(); + for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { + ((TH3F*)(fResolY->At(i)))->Reset(); + ((TH3F*)(fResolZ->At(i)))->Reset(); + ((TH3F*)(fRMSY->At(i)))->Reset(); + ((TH3F*)(fRMSZ->At(i)))->Reset(); + } + + // merge data members + cout << "histogram's merge-functins are called... " << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1 + fDeltaY->Merge(deltaYList); + fDeltaZ->Merge(deltaZList); + + TObjArray* objarray = 0; + TH1* hist = 0; + TList* histList = 0; + TIterator *objListIterator = 0; + + cout << "merging fArrayAmpRows..." << endl; + // merge fArrayAmpRows + for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) { // loop over data member, i<72 + objListIterator = arrayAmpRowList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile + if (!objarray) continue; + hist = (TProfile*)(objarray->At(i)); + histList->Add(hist); + } + ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + cout << "merging fArrayAmps..." << endl; + // merge fArrayAmps + for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72 + TIterator *objListIterator = arrayAmpList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F + if (!objarray) continue; + hist = (TH1F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH1F*)(fArrayAmp->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + cout << "merging fArrayQDY..." << endl; + // merge fArrayQDY + for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQDYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQDY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + cout << "merging fArrayQDZ..." << endl; + // merge fArrayQDZ + for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQDZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + cout << "merging fArrayQRMSY..." << endl; + // merge fArrayQRMSY + for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQRMSYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + cout << "merging fArrayQRMSZ..." << endl; + // merge fArrayQRMSZ + for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQRMSZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl; + // merge fResolY + for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = resolYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fResolY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + // merge fResolZ + for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = resolZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fResolZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + // merge fRMSY + for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = rMSYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fRMSY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + // merge fRMSZ + for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = rMSZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fRMSZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + delete deltaYList; + delete deltaZList; + delete arrayAmpRowList; + delete arrayAmpList; + delete arrayQDYList; + delete arrayQDZList; + delete arrayQRMSYList; + delete arrayQRMSZList; + delete resolYList; + delete resolZList; + delete rMSYList; + delete rMSZList; + delete nRowsList; + delete nSectList; + delete fileNoList; + delete listIterator; + + cout << "merging done!" << endl; + + return 1; +} + + +AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){ + // + // function to test AliTPCcalibTrack::Merge: + // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks" + // this object is appended 'nCalTracks' times to a TList + // A new AliTPCcalibTrack object is created which merges the list + // this object is returned + // + /* + .L AliTPCcalibTracks.cxx+g + TFile f("Output.root"); + AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks"); + //f.Close(); + TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root"); + AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param"); + clusterParamFile.Close(); + + AliTPCcalibTracks::TestMerge(calTracks, clusterParam); + */ + TList *list = new TList(); + if (ct == 0 || clusterParam == 0) return 0; + cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl; + for (Int_t i = 0; i < nCalTracks; i++) { + list->Add(new AliTPCcalibTracks(ct)); + if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl; + } + + // only for check at the end + AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(ct); + Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries(); +// Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries(); + + cout << "The list contains " << list->GetEntries() << " entries. " << endl; + + + AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018); + AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts); + cal->Merge(list); + + cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl; + cal->GetfArrayAmpRow()->At(5)->Print(); + Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries(); + + cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl; + cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " << calEntries << endl; + printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n", + calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries)); + + return cal; + +} + + + diff --git a/TPC/TPCcalib/AliTPCcalibTracks.h b/TPC/TPCcalib/AliTPCcalibTracks.h index e938c1ea8ac..836b1de8862 100644 --- a/TPC/TPCcalib/AliTPCcalibTracks.h +++ b/TPC/TPCcalib/AliTPCcalibTracks.h @@ -3,6 +3,14 @@ #include +#include +#include +#include +#include +#include +using namespace std; + +// #include "AliTPCClusterParam.h" class AliTPCClusterParam; @@ -13,57 +21,210 @@ class AliESDtrack; class TH3F; class TH1F; class TH1I; +class TTreeSRedirector; +class AliTPCcalibTracksCuts; class AliTPCcalibTracks : public TNamed { public : - - // List of branches - AliTPCcalibTracks(); - virtual ~AliTPCcalibTracks() {;} - virtual void ProofSlaveBegin(TList * output); - virtual void ProcessTrack(AliTPCseed * seed); - // - // - // + AliTPCcalibTracks(AliTPCcalibTracks* ct); + AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts); + virtual ~AliTPCcalibTracks(); + + void Process(AliTPCseed *track, AliESDtrack *esd); + Float_t TPCBetheBloch(Float_t bg); - Bool_t AcceptTrack(AliTPCseed * track); // + Bool_t AcceptTrack(AliTPCseed * track); void FillHistoCluster(AliTPCseed * track); void FillResolutionHistoLocal(AliTPCseed * track); void AlignUpDown(AliTPCseed * track, AliESDtrack *esd); - static Int_t GetBin(Float_t q,Int_t pad); - static Int_t GetBin(Int_t iq,Int_t pad); + static Int_t GetBin(Float_t q, Int_t pad); + static Int_t GetBin(Int_t iq, Int_t pad); static Float_t GetQ(Int_t bin); static Float_t GetPad(Int_t bin){return bin%3;} + Long64_t Merge(TCollection *li); + static AliTPCcalibTracks* TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks = 50); + void MakeReport(Int_t stat, char* pathName = "plots"); + void MakeAmpPlots(Int_t stat, char* pathName = "plots"); + void MakeDeltaPlots(char* pathName = "plots"); + void FitResolutionNew(char* pathName = "plots"); + void FitRMSNew(char* pathName = "plots"); + void MakeResPlotsQ(Int_t minEntries = 1, Bool_t bDraw=kFALSE, char* pathName = "plots"); + void MakeResPlotsQTree(Int_t minEntries = 1, char* pathName = "plots"); + void Draw(Option_t* opt); + void SetStyle(); + + static TH2D *MakeDiff(TH2D * hfit, TF2 * func); + static TObjArray *FitProjections(TH3F * hfit, Int_t val = 2, Int_t minEntry=500, Bool_t bDraw=kFALSE); + +//protected: + TObjArray* GetfArrayAmpRow() {return fArrayAmpRow;} + TObjArray* GetfArrayAmp() {return fArrayAmp;} + TObjArray* GetfArrayQDY() {return fArrayQDY;} + TObjArray* GetfArrayQDZ() {return fArrayQDZ;} + TObjArray* GetfArrayQRMSY() {return fArrayQRMSY;} + TObjArray* GetfArrayQRMSZ() {return fArrayQRMSZ;} + TH1F* GetfDeltaY() {return fDeltaY;} + TH1F* GetfDeltaZ() {return fDeltaZ;} + TObjArray* GetfResolY() {return fResolY;} + TObjArray* GetfResolZ() {return fResolZ;} + TObjArray* GetfRMSY() {return fRMSY;} + TObjArray* GetfRMSZ() {return fRMSZ;} + TH1I* GetfHclus() {return fHclus;} + AliTPCcalibTracksCuts* GetCuts() {return fCuts;} + + + private: - AliTPCClusterParam *fClusterParam; //pointer to cluster parameterization - TTreeSRedirector *fDebugStream; //debug stream for - TList *fOutput; //output list - // - TObjArray * fArrayAmpRow;//array with amplitudes versus row for given sector - TObjArray * fArrayAmp; //array with amplitude for sectors - TObjArray * fArrayQDY; //q binned delta Y histograms - TObjArray * fArrayQDZ; //q binned delta Z histograms - TObjArray * fArrayQRMSY; //q binned delta Y histograms - TObjArray * fArrayQRMSZ; //q binned delta Z histograms - TH1F * fDeltaY; // integrated delta y histo - TH1F * fDeltaZ; // integrated delta z histo - TObjArray * fResolY; // array of resolution histograms Y - TObjArray * fResolZ; // array of resolution histograms Z - TObjArray * fRMSY; // array of RMS histograms Y - TObjArray * fRMSZ; // array of RMS histograms Z - // - TH1I *fHclus; //! + AliTPCClusterParam *fClusterParam; //! pointer to cluster parameterization + TTreeSRedirector *fDebugStream; //! debug stream for AliTPCROC *fROC; //! - Int_t fNRows; //! - Int_t fNSect; //! - Int_t fFileNo; //! - - ClassDef(AliTPCcalibTracks,1); + TObjArray *fArrayAmpRow; // array with amplitudes versus row for given sector + TObjArray *fArrayAmp; // array with amplitude for sectors + TObjArray *fArrayQDY; // q binned delta Y histograms + TObjArray *fArrayQDZ; // q binned delta Z histograms + TObjArray *fArrayQRMSY; // q binned delta Y histograms + TObjArray *fArrayQRMSZ; // q binned delta Z histograms + TH1F *fDeltaY; // integrated delta y histo + TH1F *fDeltaZ; // integrated delta z histo + TObjArray *fResolY; // array of resolution histograms Y + TObjArray *fResolZ; // array of resolution histograms Z + TObjArray *fRMSY; // array of RMS histograms Y + TObjArray *fRMSZ; // array of RMS histograms Z + AliTPCcalibTracksCuts *fCuts; // object with cuts, that is passed to the constructor + TH1I *fHclus; // number of clusters per track + +protected: + ClassDef(AliTPCcalibTracks,1) }; +#endif + + + + +#ifndef AliTPCCALIBTRACKSCUTS_H +#define AliTPCCALIBTRACKSCUTS_H + +class AliTPCcalibTracksCuts: public TNamed { + ////////////////////////////////////////////////////// + // // + // Class to specify cuts for track analysis // + // with AliTPCcalibTracks // + // // + ////////////////////////////////////////////////////// + +public: + AliTPCcalibTracksCuts(Int_t minClusters = 20, Float_t minRatio = 0.4, Float_t max1pt = 0.5, + Float_t edgeXZCutNoise = 0.13, Float_t edgeThetaCutNoise = 0.018): + TNamed("calibTracksCuts", "calibTracksCuts") { + // + // Constuctor for AliTPCcalibTracksCuts + // specify the cuts to be set on the processed tracks + // + fMinClusters = minClusters; + fMinRatio = minRatio; + fMax1pt = max1pt; + fEdgeYXCutNoise = edgeXZCutNoise; + fEdgeThetaCutNoise = edgeThetaCutNoise; + } + virtual ~AliTPCcalibTracksCuts(){cout << "AliTPCcalibTracksCuts destructor called, nothing happend." << endl;} + void SetMinClusters(Int_t minClusters){fMinClusters = minClusters;} + void SetMinRatio(Float_t minRatio){fMinRatio = minRatio;} + void SetMax1pt(Float_t max1pt){fMax1pt = max1pt;} + void SetEdgeXYCutNoise(Float_t edgeCutNoise){fEdgeYXCutNoise = edgeCutNoise;} + void SetEdgeThetaCutNoise(Float_t edgeCutNoise){fEdgeThetaCutNoise = edgeCutNoise;} + Int_t GetMinClusters(){return fMinClusters;} + Float_t GetMinRatio(){return fMinRatio;} + Float_t GetMax1pt(){return fMax1pt;} + Float_t GetEdgeYXCutNoise(){return fEdgeYXCutNoise;} + Float_t GetEdgeThetaCutNoise(){return fEdgeThetaCutNoise;} + void Print(){ + cout << ": The following cuts are specified: " << endl; + cout << "fMinClusters: " << fMinClusters << endl; + cout << "fMinRatio: " << fMinRatio << endl; + cout << "fMax1pt: " << fMax1pt << endl; + cout << "fEdgeYXCutNoise: " << fEdgeYXCutNoise << endl; + cout << "fEdgeThetaCutNoise: " << fEdgeThetaCutNoise << endl; + } // Prints out the specified cuts + +private: + Int_t fMinClusters; // number of clusters + Float_t fMinRatio; // kMinRratio = 0.4 + Float_t fMax1pt; // kMax1pt = 0.5 + Float_t fEdgeYXCutNoise; // kEdgeYXCutNoise = 0.13 + Float_t fEdgeThetaCutNoise; // kEdgeThetaCutNoise = 0.018 + +protected: + ClassDef(AliTPCcalibTracksCuts,1) +}; #endif + + +#ifndef AliTPCCALIBTRACKSCUTS_H +#define AliTPCCALIBTRACKSCUTS_H + + +class AliTPCcalibTracksCuts: public TNamed { + ////////////////////////////////////////////////////// + // // + // Class to specify cuts for track analysis // + // with AliTPCcalibTracks // + // // + ////////////////////////////////////////////////////// + +public: + AliTPCcalibTracksCuts(Int_t minClusters = 20, Float_t minRatio = 0.4, Float_t max1pt = 0.5, + Float_t edgeXZCutNoise = 0.13, Float_t edgeThetaCutNoise = 0.018): + TNamed("calibTracksCuts", "calibTracksCuts") { + // + // Constuctor for AliTPCcalibTracksCuts + // specify the cuts to be set on the processed tracks + // + fMinClusters = minClusters; + fMinRatio = minRatio; + fMax1pt = max1pt; + fEdgeYXCutNoise = edgeXZCutNoise; + fEdgeThetaCutNoise = edgeThetaCutNoise; + } + virtual ~AliTPCcalibTracksCuts(){cout << "AliTPCcalibTracksCuts destructor called, nothing happend." << endl;} + void SetMinClusters(Int_t minClusters){fMinClusters = minClusters;} + void SetMinRatio(Float_t minRatio){fMinRatio = minRatio;} + void SetMax1pt(Float_t max1pt){fMax1pt = max1pt;} + void SetEdgeXYCutNoise(Float_t edgeCutNoise){fEdgeYXCutNoise = edgeCutNoise;} + void SetEdgeThetaCutNoise(Float_t edgeCutNoise){fEdgeThetaCutNoise = edgeCutNoise;} + Int_t GetMinClusters(){return fMinClusters;} + Float_t GetMinRatio(){return fMinRatio;} + Float_t GetMax1pt(){return fMax1pt;} + Float_t GetEdgeYXCutNoise(){return fEdgeYXCutNoise;} + Float_t GetEdgeThetaCutNoise(){return fEdgeThetaCutNoise;} +// void Print(); // Prints out the specified cuts + void Print(){ + cout << ": The following cuts are specified: " << endl; + cout << "fMinClusters: " << fMinClusters << endl; + cout << "fMinRatio: " << fMinRatio << endl; + cout << "fMax1pt: " << fMax1pt << endl; + cout << "fEdgeYXCutNoise: " << fEdgeYXCutNoise << endl; + cout << "fEdgeThetaCutNoise: " << fEdgeThetaCutNoise << endl; + } // Prints out the specified cuts + + +private: + Int_t fMinClusters; // number of clusters + Float_t fMinRatio; // kMinRratio = 0.4 + Float_t fMax1pt; // kMax1pt = 0.5 + Float_t fEdgeYXCutNoise; // kEdgeYXCutNoise = 0.13 + Float_t fEdgeThetaCutNoise; // kEdgeThetaCutNoise = 0.018 + +protected: + ClassDef(AliTPCcalibTracksCuts,1) +}; + + +#endif + + diff --git a/TPC/TPCcalib/AliTPCcalibV0.cxx b/TPC/TPCcalib/AliTPCcalibV0.cxx index 867bca37795..bcdf3613807 100644 --- a/TPC/TPCcalib/AliTPCcalibV0.cxx +++ b/TPC/TPCcalib/AliTPCcalibV0.cxx @@ -63,6 +63,16 @@ AliTPCcalibV0::AliTPCcalibV0() : G__SetCatchException(0); fDebugStream = new TTreeSRedirector("V0debug.root"); fPdg = new TDatabasePDG; + + + // create output histograms + fTPCdEdx = new TH2F("TPCdEdX", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + BinLogX(fTPCdEdx); + fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + fTPCdEdxP = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + + } AliTPCcalibV0::~AliTPCcalibV0(){ @@ -155,14 +165,18 @@ void AliTPCcalibV0::MakeMC(){ if (p0->Pt()Vz()>250) findable= kFALSE; if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; - if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++; - + if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE; + else + if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++; + p1 = fStack->Particle(id1); if (p1->R()>kMaxRad) findable = kFALSE; if (p1->Pt()Vz())>250) findable= kFALSE; if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; - if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++; + if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE; + else + if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++; } // (*fDebugStream)<<"MC0"<< @@ -220,8 +234,8 @@ void AliTPCcalibV0::MakeMC(){ "Pp.="<AddLast(v0); +// // +// // +// // +// delete v0K0; +// delete v0Gamma; +// delete v0Lambda42; +// delete v0Lambda24; +// delete v0K0C; +// delete v0GammaC; +// delete v0Lambda42C; +// delete v0Lambda24C; +// } +// ProcessPI0(); +// } + + + + void AliTPCcalibV0::ProcessV0(Int_t ftype){ // // - const Double_t ktimeK0 = 2.684; - const Double_t ktimeLambda = 7.89; - if (! fGammas) fGammas = new TObjArray(10); fGammas->Clear(); @@ -333,15 +505,21 @@ void AliTPCcalibV0::ProcessV0(Int_t ftype){ // for (Int_t ivertex=0; ivertexAt(ivertex); - AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); - AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); + AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track + AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track + + const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam(); + const AliExternalTrackParam * paramInPos = trackP->GetInnerParam(); + + if (!paramInPos) continue; // in case the inner paramters do not exist + if (!paramInNeg) continue; // // // - AliKFParticle *v0K0 = Fit(primVtx,v0,211,211); + AliKFParticle *v0K0 = Fit(primVtx,v0,-211,211); AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11); - AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211); - AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212); + AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211); + AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212); //Set production vertex v0K0->SetProductionVertex( primVtx ); v0Gamma->SetProductionVertex( primVtx ); @@ -359,10 +537,10 @@ void AliTPCcalibV0::ProcessV0(Int_t ftype){ // // Mass Contrained params // - AliKFParticle *v0K0C = Fit(primVtx,v0,211,211); + AliKFParticle *v0K0C = Fit(primVtx,v0,-211,211); AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); - AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211); - AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212); + AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar + AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda // v0K0C->SetProductionVertex( primVtx ); v0GammaC->SetProductionVertex( primVtx ); @@ -371,8 +549,8 @@ void AliTPCcalibV0::ProcessV0(Int_t ftype){ v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); v0GammaC->SetMassConstraint(0); - v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); - v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); + v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); + v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); // Double_t timeK0, sigmaTimeK0; Double_t timeLambda42, sigmaTimeLambda42; @@ -405,14 +583,49 @@ void AliTPCcalibV0::ProcessV0(Int_t ftype){ if (chi2GammaGetP()/fPdg->GetParticle(-2212)->Mass(); - Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass(); - Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass(); - Float_t dedxTeorP = TPCBetheBloch(betaGammaP); - Float_t dedxTeorPi = TPCBetheBloch(betaGammaPi);; - Float_t dedxTeorEl = TPCBetheBloch(betaGammaEl);; + + // 0 is negative particle; 1 is positive particle + Float_t betaGamma0 = 0; + Float_t betaGamma1 = 0; + + switch (type) { + case 0: + betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); + break; + case 1: + betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass(); + break; + case 2: + betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); + break; + case 3: + betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass(); + break; + } + + // cuts and histogram filling + Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2 + + if (minChi2C < 2 && ftype == 1) { + // + if (chi2K0C < 10*minChi2C) numCand++; + if (chi2GammaC < 10*minChi2C) numCand++; + if (chi2Lambda42C < 10*minChi2C) numCand++; + if (chi2Lambda24C < 10*minChi2C) numCand++; + // + if (numCand < 2) { + if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal()); + if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal()); + } + } + // // + // write output tree if (minChi2>50) continue; (*fDebugStream)<<"V0"<< "ftype="<