X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCcalibTracks.cxx;h=39d8f4919eec5828bbbbecc6dd1b025c0ea859a0;hb=a054a582ed435c7bbc2680fe557b57c32ba19333;hp=08d3e431c10b282d27fb9d074f4c40a11ec2e79e;hpb=2c632057b20a755472de41d084d02d410a4d451d;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCcalibTracks.cxx b/TPC/AliTPCcalibTracks.cxx index 08d3e431c10..39d8f4919ee 100644 --- a/TPC/AliTPCcalibTracks.cxx +++ b/TPC/AliTPCcalibTracks.cxx @@ -16,12 +16,49 @@ /////////////////////////////////////////////////////////////////////////////// // // -// 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. // +// Class for calibration of the cluster properties: // +// Cluster position resolution (RMS) and short term distortions (within pad or within time bin) +// +// Algorithm: +// 1. Creation of the residual/properties histograms +// 2. Filling of the histograms. +// 2.a Traklet fitting +// 2.b Resudual filling +// 3. Postprocessing: Creation of the tree with the mean values and RMS at different bins +// 4. : Paramaterization fitting. +// In old AliRoot the RMS paramterization was used to characterize cluster resolution +// Mean value /bias was neglected +// 5. To be implemented +// a.) lookup table for the distortion parmaterization - THn +// b.) Usage in the reconstruction +// +// +// 1. Creation of the histograms: MakeHistos() +// +// 6 n dimensional histograms are filled during the calibration: +// cluster performance bins +// 0 - variable of interest +// 1 - pad type - 0- IROC Short 1- OCROC medium 2 - OROC long pads +// 2 - drift length - drift length -0-1 +// 3 - Qmax - Qmax - 2- 400 +// 4 - cog - COG position - 0-1 +// 5 - tan(phi) - local angle +// Histograms: +// THnSparse *fHisDeltaY; // THnSparse - delta Y between the cluster and tracklet interpolation (+-X(5?) rows) +// THnSparse *fHisDeltaZ; // THnSparse - delta Z +// THnSparse *fHisRMSY; // THnSparse - rms Y +// THnSparse *fHisRMSZ; // THnSparse - rms Z +// THnSparse *fHisQmax; // THnSparse - qmax +// THnSparse *fHisQtot; // THnSparse - qtot +// +// +// + + + // // 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. // @@ -34,6 +71,24 @@ Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam) */ +/* + Example usage - creation of the residual trees: + TFile f("CalibObjects.root"); + AliTPCcalibTracks *calibTracks = ( AliTPCcalibTracks *)TPCCalib->FindObject("calibTracks"); + TH2 * his2 = calibTracks->fHisDeltaY->Projection(0,4); + his2->SetName("his2"); + his2->FitSlicesY(); + + + TTreeSRedirector *pcstream = new TTreeSRedirector("clusterLookup.root"); + AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaY, pcstream, 0); + AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaZ, pcstream, 1); + delete pcstream; + TFile fl("clusterLookup.root"); + TCut cutNI="cogA==0&&angleA==0&&qmaxA==0&&zA==0&&ipadA==0" + +*/ + // /////////////////////////////////////////////////////////////////////////////// @@ -48,6 +103,8 @@ using namespace std; #include #include #include +#include + // //#include #include @@ -66,6 +123,7 @@ using namespace std; #include #include #include +#include // // AliROOT includes @@ -85,201 +143,168 @@ using namespace std; #include "AliTPCcalibTracks.h" #include "AliTPCClusterParam.h" #include "AliTPCcalibTracksCuts.h" -#include "AliTPCCalPadRegion.h" #include "AliTPCCalPad.h" #include "AliTPCCalROC.h" #include "TText.h" #include "TPaveText.h" #include "TSystem.h" +#include "TStatToolkit.h" +#include "TCut.h" +#include "THnSparse.h" +#include "AliRieman.h" +#include "TRandom.h" -// Thread-stuff -//#include "TThread.h" -//#include "TMutex.h" -//#include "TLockFile.h" +Double_t AliTPCcalibTracks::fgkMergeEntriesCut=10000000.; //10**7 clusters ClassImp(AliTPCcalibTracks) AliTPCcalibTracks::AliTPCcalibTracks(): - TNamed(), + AliTPCcalibBase(), fClusterParam(0), - fDebugStream(0), fROC(0), - fArrayAmpRow(0), - fArrayAmp(0), + fHisDeltaY(0), // THnSparse - delta Y + fHisDeltaZ(0), // THnSparse - delta Z + fHisRMSY(0), // THnSparse - rms Y + fHisRMSZ(0), // THnSparse - rms Z + fHisQmax(0), // THnSparse - qmax + fHisQtot(0), // THnSparse - qtot + fPtDownscaleRatio(20), // pt downscaling ratio (use subsample of data) + fQDownscaleRatio(20), // Q downscaling ratio (use subsample of dta) fArrayQDY(0), fArrayQDZ(0), fArrayQRMSY(0), fArrayQRMSZ(0), - fArrayChargeVsDriftlength(0), - fcalPadRegionChargeVsDriftlength(0), - fDeltaY(0), - fDeltaZ(0), fResolY(0), fResolZ(0), fRMSY(0), fRMSZ(0), fCuts(0), - fHclus(0), fRejectedTracksHisto(0), - fHclusterPerPadrow(0), - fHclusterPerPadrowRaw(0), fClusterCutHisto(0), fCalPadClusterPerPad(0), - fCalPadClusterPerPadRaw(0), - fDebugLevel(0), - fFitterLinY1(0), //! - fFitterLinZ1(0), //! - fFitterLinY2(0), //! - fFitterLinZ2(0), //! - fFitterParY(0), //! - fFitterParZ(0) //! + fCalPadClusterPerPadRaw(0) { // // AliTPCcalibTracks default constructor // - if (fDebugLevel > 0) cout << "AliTPCcalibTracks' default constructor called" << endl; + if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl; } -AliTPCcalibTracks::AliTPCcalibTracks(AliTPCcalibTracks* ct): -TNamed(), + +AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): + AliTPCcalibBase(calibTracks), fClusterParam(0), - fDebugStream(0), fROC(0), - fArrayAmpRow(0), - fArrayAmp(0), + fHisDeltaY(0), // THnSparse - delta Y + fHisDeltaZ(0), // THnSparse - delta Z + fHisRMSY(0), // THnSparse - rms Y + fHisRMSZ(0), // THnSparse - rms Z + fHisQmax(0), // THnSparse - qmax + fHisQtot(0), // THnSparse - qtot + fPtDownscaleRatio(20), // pt downscaling ratio (use subsample of data) + fQDownscaleRatio(20), // Q downscaling ratio (use subsample of dta) fArrayQDY(0), fArrayQDZ(0), fArrayQRMSY(0), fArrayQRMSZ(0), - fArrayChargeVsDriftlength(0), - fcalPadRegionChargeVsDriftlength(0), - fDeltaY(0), - fDeltaZ(0), fResolY(0), fResolZ(0), fRMSY(0), fRMSZ(0), fCuts(0), - fHclus(0), fRejectedTracksHisto(0), - fHclusterPerPadrow(0), - fHclusterPerPadrowRaw(0), fClusterCutHisto(0), fCalPadClusterPerPad(0), - fCalPadClusterPerPadRaw(0), - fDebugLevel(0), - fFitterLinY1(0), //! - fFitterLinZ1(0), //! - fFitterLinY2(0), //! - fFitterLinZ2(0), //! - fFitterParY(0), //! - fFitterParZ(0) //! + fCalPadClusterPerPadRaw(0) { // // AliTPCcalibTracks copy constructor // - if (fDebugLevel > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl; + if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl; Bool_t dirStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); Int_t length = -1; - // backward compatibility: if the data member doesn't yet exist, it will not be merged - (ct->fArrayAmpRow) ? length = ct->fArrayAmpRow->GetEntriesFast() : length = -1; - 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); - } - (ct->fArrayQDY) ? length = ct->fArrayQDY->GetEntriesFast() : length = -1; + (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1; 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); + fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i); + fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i); + fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i); + fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i); } - (ct->fResolY) ? length = ct->fResolY->GetEntriesFast() : length = -1; + (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1; 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); + fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i); + fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i); + fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i); + fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i); } - (ct->fArrayChargeVsDriftlength) ? length = ct->fArrayChargeVsDriftlength->GetEntriesFast() : length = -1; - (ct->fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0; - for (Int_t i = 0; i < length; i++) { - fArrayChargeVsDriftlength->AddAt( ((TProfile*)ct->fArrayChargeVsDriftlength->At(i)->Clone()), i); - } - fDeltaY = (TH1F*)ct->fDeltaY->Clone(); - fDeltaZ = (TH1F*)ct->fDeltaZ->Clone(); - fHclus = (TH1I*)ct->fHclus->Clone(); - fClusterCutHisto = (TH2I*)ct->fClusterCutHisto->Clone(); - fRejectedTracksHisto = (TH1I*)ct->fRejectedTracksHisto->Clone(); - fHclusterPerPadrow = (TH1I*)ct->fHclusterPerPadrow->Clone(); - fHclusterPerPadrowRaw = (TH1I*)ct->fHclusterPerPadrowRaw->Clone(); - fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)ct->fcalPadRegionChargeVsDriftlength->Clone(); - fCalPadClusterPerPad = (AliTPCCalPad*)ct->fCalPadClusterPerPad->Clone(); - fCalPadClusterPerPadRaw = (AliTPCCalPad*)ct->fCalPadClusterPerPadRaw->Clone(); - - fCuts = new AliTPCcalibTracksCuts(ct->fCuts->GetMinClusters(), ct->fCuts->GetMinRatio(), - ct->fCuts->GetMax1pt(), ct->fCuts->GetEdgeYXCutNoise(), ct->fCuts->GetEdgeThetaCutNoise()); - fDebugLevel = ct->GetLogLevel(); - SetNameTitle(ct->GetName(), ct->GetTitle()); + fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone(); + fRejectedTracksHisto = (TH1I*)calibTracks.fRejectedTracksHisto->Clone(); + fCalPadClusterPerPad = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone(); + fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone(); + + fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(), + calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise()); + SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle()); TH1::AddDirectory(dirStatus); // set status back to original status // cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED } +AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){ + // + // assgnment operator + // + if (this != &calibTracks) { + new (this) AliTPCcalibTracks(calibTracks); + } + return *this; + +} + + AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) : - TNamed(name, title), + AliTPCcalibBase(), fClusterParam(0), - fDebugStream(0), fROC(0), - fArrayAmpRow(0), - fArrayAmp(0), + fHisDeltaY(0), // THnSparse - delta Y + fHisDeltaZ(0), // THnSparse - delta Z + fHisRMSY(0), // THnSparse - rms Y + fHisRMSZ(0), // THnSparse - rms Z + fHisQmax(0), // THnSparse - qmax + fHisQtot(0), // THnSparse - qtot + fPtDownscaleRatio(20), // pt downscaling ratio (use subsample of data) + fQDownscaleRatio(20), // Q downscaling ratio (use subsample of dta) fArrayQDY(0), fArrayQDZ(0), fArrayQRMSY(0), fArrayQRMSZ(0), - fArrayChargeVsDriftlength(0), - fcalPadRegionChargeVsDriftlength(0), - fDeltaY(0), - fDeltaZ(0), fResolY(0), fResolZ(0), fRMSY(0), fRMSZ(0), fCuts(0), - fHclus(0), fRejectedTracksHisto(0), - fHclusterPerPadrow(0), - fHclusterPerPadrowRaw(0), fClusterCutHisto(0), fCalPadClusterPerPad(0), - fCalPadClusterPerPadRaw(0), - fDebugLevel(0), - fFitterLinY1(0), //! - fFitterLinZ1(0), //! - fFitterLinY2(0), //! - fFitterLinZ2(0), //! - fFitterParY(0), //! - fFitterParZ(0) //! + fCalPadClusterPerPadRaw(0) { // // AliTPCcalibTracks constructor @@ -287,11 +312,15 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al // 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. - // log level - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStream, 6: waste your screen + // + // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen // // All histograms are instatiated in this constructor. // - if (fDebugLevel > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl; + this->SetName(name); + this->SetTitle(title); + + if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl; G__SetCatchException(0); fClusterParam = clusterParam; @@ -302,66 +331,19 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)"); } fCuts = cuts; - fDebugLevel = logLevel; - if (fDebugLevel > 4) fDebugStream = new TTreeSRedirector("TPCSelectorDebug.root"); // needs investigation !!!!! + SetDebugLevel(logLevel); + MakeHistos(); TH1::AddDirectory(kFALSE); - char chname[1000]; - TProfile * prof1=0; - TH1F * his1 =0; - fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160); // valgrind 3 - fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 10, 1, 10); - fHclusterPerPadrow = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160); - fHclusterPerPadrowRaw = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160); + fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10); fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad"); fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters"); fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159); - // Amplitude - sector - row histograms - fArrayAmpRow = new TObjArray(72); - fArrayAmp = new TObjArray(72); - fArrayChargeVsDriftlength = 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); - - // driftlength - sprintf(chname, "driftlengt vs. charge, ROC %i", i); - prof1 = new TProfile(chname, chname, 500, 0, 250); - prof1->SetYTitle("Charge"); - prof1->SetXTitle("Driftlength"); - fArrayChargeVsDriftlength->AddAt(prof1,i); - sprintf(chname, "driftlengt vs. charge, ROC %i", i+36); - prof1 = new TProfile(chname, chname, 500, 0, 250); - prof1->SetYTitle("Charge"); - prof1->SetXTitle("Driftlength"); - fArrayChargeVsDriftlength->AddAt(prof1,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); @@ -408,44 +390,25 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al 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); + char hname[200]; + snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -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); + snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -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); + snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6); 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); + snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6); fArrayQRMSZ->AddAt(his3D, bin); } } - fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region"); - TProfile *tempProf; - for (UInt_t padSize = 0; padSize < 3; padSize++) { - for (UInt_t isector = 0; isector < 36; isector++) { - if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector); - if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium pads", isector); - if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long pads", isector); - tempProf = new TProfile(chname, chname, 500, 0, 250); - tempProf->SetYTitle("Charge"); - tempProf->SetXTitle("Driftlength"); - fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize); - } - } - fFitterLinY1 = new TLinearFitter (2,"pol1"); - fFitterLinZ1 = new TLinearFitter (2,"pol1"); - fFitterLinY2 = new TLinearFitter (2,"pol1"); - fFitterLinZ2 = new TLinearFitter (2,"pol1"); - fFitterParY = new TLinearFitter (3,"pol2"); - fFitterParZ = new TLinearFitter (3,"pol2"); - - if (fDebugLevel > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; + + if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; cout << "end of main constructor" << endl; // TO BE REMOVED } @@ -455,94 +418,71 @@ AliTPCcalibTracks::~AliTPCcalibTracks() { // AliTPCcalibTracks destructor // - if (fDebugLevel > 0) cout << "AliTPCcalibTracks' destuctor called." << endl; + if (GetDebugLevel() > 0) 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); + 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 (fArrayChargeVsDriftlength) length = fArrayChargeVsDriftlength->GetEntriesFast(); - for (Int_t i = 0; i < length; i++){ - delete fArrayChargeVsDriftlength->At(i); + 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 fFitterLinY1; - delete fFitterLinZ1; - delete fFitterLinY2; - delete fFitterLinZ2; - delete fFitterParY; - delete fFitterParZ; + delete fArrayQDY; delete fArrayQDZ; delete fArrayQRMSY; delete fArrayQRMSZ; - delete fArrayChargeVsDriftlength; - delete fHclus; delete fRejectedTracksHisto; delete fClusterCutHisto; - delete fHclusterPerPadrow; - delete fHclusterPerPadrowRaw; if (fCalPadClusterPerPad) delete fCalPadClusterPerPad; if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw; - fcalPadRegionChargeVsDriftlength->Delete(); - delete fcalPadRegionChargeVsDriftlength; - if (fDebugLevel > 4) delete fDebugStream; + delete fHisDeltaY; // THnSparse - delta Y + delete fHisDeltaZ; // THnSparse - delta Z + delete fHisRMSY; // THnSparse - rms Y + delete fHisRMSZ; // THnSparse - rms Z + delete fHisQmax; // THnSparse - qmax + delete fHisQtot; // THnSparse - qtot + } -void AliTPCcalibTracks::AddInfo(TChain * chain, char* fileName){ - // - // Add the neccessary information for processing to the chain - // (cluster parametrization) - // - TFile clusterParamFile(fileName); - AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param"); - chain->GetUserInfo()->AddLast((TObject*)clusterParam); - cout << "Clusterparametrization added to the chain." << endl; -} - -void AliTPCcalibTracks::Process(AliTPCseed *track, AliESDtrack *esd){ +void AliTPCcalibTracks::Process(AliTPCseed *track){ // // To be called in the selector // first AcceptTrack is evaluated, then calls all the following analyse functions: // FillResolutionHistoLocal(track) - // AlignUpDown(track, esd) + // - if (fDebugLevel > 5) Info("Process","Starting to process the track..."); + Double_t scalept= TMath::Min(1./TMath::Abs(track->GetParameter()[4]),2.)/0.5; + Bool_t isSelected = (TMath::Exp(scalept)>fPtDownscaleRatio*gRandom->Rndm()); + if (!isSelected) return; + + if (GetDebugLevel() > 5) Info("Process","Starting to process the track..."); Int_t accpetStatus = AcceptTrack(track); if (accpetStatus == 0) { FillResolutionHistoLocal(track); - // AlignUpDown(track, esd); } else fRejectedTracksHisto->Fill(accpetStatus); } @@ -619,7 +559,7 @@ Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){ //if (TMath::Abs(track->GetZ())<10.) return kFALSE; //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE; - if (fDebugLevel > 5) Info("AcceptTrack","Track has been accepted."); + if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted."); return 0; } @@ -642,616 +582,328 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ // 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 + // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY // // write debug information to 'TPCSelectorDebug.root' // only for every kDeltaWriteDebugStream'th padrow to reduce data volume // and to avoid redundant data // - if (fDebugLevel > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****"); - 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 - const Int_t kDeltaWriteDebugStream = 5; // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream -// TLinearFitter fFitterLinY1 = fFitterLinY1; -// TLinearFitter fFitterLinZ1 = ffFitterLinZ1; -// TLinearFitter fFitterLinY2 = ffFitterLinY2; -// TLinearFitter fFitterLinZ2 = ffFitterLinZ2; -// TLinearFitter fFitterParY = ffFitterParY; -// TLinearFitter fFitterParZ = ffFitterParZ; - 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; // number of tracklets for given track - Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction - Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction - Int_t nClusters = 0; // working variable, number of clusters per tracklet - Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet - - fHclus->Fill(track->GetNumberOfClusters()); // for statistics overview - // --------------------------------------------------------------------- - 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; // no cluster found - Int_t sector = cluster0->GetDetector(); - fHclusterPerPadrowRaw->Fill(irow); - - Int_t ipad = TMath::Nint(cluster0->GetPad()); - Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad())); - fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 ); - - if (sector != sectorG){ - // track leaves sector before it crossed enough rows to fit / initialization - nClusters = 0; - fFitterParY->ClearPoints(); - fFitterParZ->ClearPoints(); - sectorG = sector; + /* {//SG + static long n=0; + if( n==0 ){ + if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){ + cout<<"Map found "<GetX(); - fFitterParY->AddPoint(&x, cluster0->GetY(), 1); - fFitterParZ->AddPoint(&x, cluster0->GetZ(), 1); - // - if ( nClusters >= kDelta + 3 ){ - // if more than 13 (kDelta+3) clusters were added to the fitters - // fit the tracklet, increase trackletCounter - fFitterParY->Eval(); - fFitterParZ->Eval(); - nTrackletsAll++; - csigmaY += fFitterParY->GetChisquare() / (nClusters - 3.); - csigmaZ += fFitterParZ->GetChisquare() / (nClusters - 3.); - nClusters = -1; - fFitterParY->ClearPoints(); - fFitterParZ->ClearPoints(); - } + } + if( n%100 ==0 )cout< 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****"); + const Int_t kDelta = 10; // delta rows to fit + const Float_t kMinRatio = 0.75; // minimal ratio + const Float_t kChi2Cut = 10.; // cut chi2 - left right + const Float_t kSigmaCut = 3.; //sigma cut + const Float_t kdEdxCut=300; + const Float_t kPtCut=0.300; + + if (track->GetTPCsignal()>kdEdxCut) return; + if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())GetClusterPointer(irow); + if (!cluster0) continue; // no cluster found + Int_t sector = cluster0->GetDetector(); + + Int_t ipad = TMath::Nint(cluster0->GetPad()); + Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad())); + fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 ); + + if (sector != sectorG){ + // track leaves sector before it crossed enough rows to fit / initialization + nClusters = 0; + fFitterParY.ClearPoints(); + fFitterParZ.ClearPoints(); + sectorG = sector; + } + else { + nClusters++; + if (refX<1) refX=cluster0->GetX()+kDelta*0.5; + Double_t x = cluster0->GetX()-refX; + fFitterParY.AddPoint(&x, cluster0->GetY(), 1); + fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1); + // + if ( nClusters >= kDelta + 3 ){ + // if more than 13 (kDelta+3) clusters were added to the fitters + // fit the tracklet, increase trackletCounter + fFitterParY.Eval(); + fFitterParZ.Eval(); + nTrackletsAll++; + csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.); + csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.); + nClusters = -1; + fFitterParY.ClearPoints(); + fFitterParZ.ClearPoints(); + refX=0; } - } // 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++){ - // loop again over all rows along the track - // do analysis - // - Int_t nclFound = 0; // number of clusters in the neighborhood - Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster - Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster - AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); - if (!cluster0) continue; - Int_t sector = cluster0->GetDetector(); - Float_t xref = cluster0->GetX(); - - // Make Fit - fFitterParY->ClearPoints(); - fFitterParZ->ClearPoints(); - fFitterLinY1->ClearPoints(); - fFitterLinZ1->ClearPoints(); - fFitterLinY2->ClearPoints(); - fFitterLinZ2->ClearPoints(); - - // 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 irow = 0; irow < 159; irow++) + // mean chi^2 for all tracklet fits in Y and in Z direction: + csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1)); + csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1)); + if (csigmaY<=TMath::KUncertainty() || csigmaZ<= TMath::KUncertainty()) return; + // --------------------------------------------------------------------- + // + // + + for (Int_t irow = kDelta; irow < 159-kDelta; irow++){ + // loop again over all rows along the track + // do analysis + // + Int_t nclFound = 0; // number of clusters in the neighborhood + Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster + Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster + AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); + if (!cluster0) continue; + Int_t sector = cluster0->GetDetector(); + Float_t xref = cluster0->GetX(); + + // Make Fit + fFitterParY.ClearPoints(); + fFitterParZ.ClearPoints(); + fFitterParYWeight.ClearPoints(); + fFitterParZWeight.ClearPoints(); + // 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 + + // first fit the track angle for wave correction + + AliRieman riemanFitAngle( 2*kDelta ); //SG + + if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){ for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ - // loop over irow +- kDelta rows (neighboured rows) - // - // - if (idelta == 0) continue; // don't use center cluster - if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC - AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta); - if (!currentCluster) continue; - if (currentCluster->GetType() < 0) continue; - if (currentCluster->GetDetector() != sector) continue; - Double_t x = currentCluster->GetX() - xref; // x = differece: current cluster - cluster @ irow - nclFound++; - if (idelta < 0){ - ncl0++; - fFitterLinY1->AddPoint(&x, currentCluster->GetY(), csigmaY); - fFitterLinZ1->AddPoint(&x, currentCluster->GetZ(), csigmaZ); - } - if (idelta > 0){ - ncl1++; - fFitterLinY2->AddPoint(&x, currentCluster->GetY(), csigmaY); - fFitterLinZ2->AddPoint(&x, currentCluster->GetZ(), csigmaZ); - } - fFitterParY->AddPoint(&x, currentCluster->GetY(), csigmaY); - fFitterParZ->AddPoint(&x, currentCluster->GetZ(), csigmaZ); - } // loop over neighbourhood for fitter filling - - if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10); - if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow); - if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow - fFitterParY->Eval(); - fFitterParZ->Eval(); - Double_t chi2 = (fFitterParY->GetChisquare() + fFitterParZ->GetChisquare()) / (2. * nclFound - 6.); - if (chi2 > kCutChi2) fRejectedTracksHisto->Fill(9); - if (chi2 > kCutChi2) fClusterCutHisto->Fill(2, irow); - if (chi2 > kCutChi2) continue; // if chi^2 is too big goto next padrow - - // REMOVE KINK - // only when there are enough clusters (4) in each direction - if (ncl0 > 4){ - fFitterLinY1->Eval(); - fFitterLinZ1->Eval(); - } - if (ncl1 > 4){ - fFitterLinY2->Eval(); - fFitterLinZ2->Eval(); + // loop over irow +- kDelta rows (neighboured rows) + if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC + AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta); + if (!currentCluster) continue; + if (currentCluster->GetType() < 0) continue; + if (currentCluster->GetDetector() != sector) continue; + riemanFitAngle.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ); } - - if (ncl0 > 4 && ncl1 > 4){ - fFitterLinY1->GetCovarianceMatrix(matrixY0); - fFitterLinY2->GetCovarianceMatrix(matrixY1); - fFitterLinZ1->GetCovarianceMatrix(matrixZ0); - fFitterLinZ2->GetCovarianceMatrix(matrixZ1); - fFitterLinY2->GetParameters(paramY1); - fFitterLinZ2->GetParameters(paramZ1); - fFitterLinY1->GetParameters(paramY0); - fFitterLinZ1->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()); - 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()); - matrixZ0.Invert(); - TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ); - TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ); - chi2 += chi2Z(0, 0); - - // REMOVE KINK - if (chi2 * 0.25 > kCutChi2) fRejectedTracksHisto->Fill(8); - if (chi2 * 0.25 > kCutChi2) fClusterCutHisto->Fill(3, irow); - 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 + if( riemanFitAngle.GetN()>3 ) riemanFitAngle.Update(); + } + + // do fit + + AliRieman riemanFit(2*kDelta); + AliRieman riemanFitW(2*kDelta); + for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ + // loop over irow +- kDelta rows (neighboured rows) + // + // + if (idelta == 0) continue; // don't use center cluster + if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC + AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta); + if (!currentCluster) continue; + if (currentCluster->GetType() < 0) continue; + if (currentCluster->GetDetector() != sector) continue; + nclFound++; + if (idelta < 0){ + ncl0++; } - - // current padrow has no kink - - // get fit parameters from pol2 fit: - Double_t paramY[4], paramZ[4]; - paramY[0] = fFitterParY->GetParameter(0); - paramY[1] = fFitterParY->GetParameter(1); - paramY[2] = fFitterParY->GetParameter(2); - paramZ[0] = fFitterParZ->GetParameter(0); - paramZ[1] = fFitterParZ->GetParameter(1); - paramZ[2] = fFitterParZ->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 anglez = paramZ[1]; - - Float_t max = cluster0->GetMax(); - UInt_t isegment = cluster0->GetDetector() % 36; - Int_t padSize = 0; // short pads - if (cluster0->GetDetector() >= 36) { - padSize = 1; // medium pads - if (cluster0->GetRow() > 63) padSize = 2; // long pads + if (idelta > 0){ + ncl1++; } - - // ========================================= - // wirte collected information to histograms - // ========================================= - - TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector); - if ( irow >= kFirstLargePad) max /= kLargePadSize; - profAmpRow->Fill( (Double_t)cluster0->GetRow(), max ); - TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector); - hisAmp->Fill(max); - - // remove the following two lines one day: - TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector); - profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max ); - - TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize)); - profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max ); - - fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow - Int_t ipad = TMath::Nint(cluster0->GetPad()); - Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad())); - fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 ); - - - TH3F * his3 = 0; - his3 = (TH3F*)fRMSY->At(padSize); - if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) ); - his3 = (TH3F*)fRMSZ->At(padSize); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) ); - - his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize)); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) ); - his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize)); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) ); - - - // Fill resolution histograms - Bool_t useForResol = kTRUE; - if (fFitterParY->GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE; - - if (useForResol){ - fDeltaY->Fill(deltay); - fDeltaZ->Fill(deltaz); - his3 = (TH3F*)fResolY->At(padSize); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay ); - his3 = (TH3F*)fResolZ->At(padSize); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz ); - his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize)); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay ); - his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize)); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz ); + //SG!!! + double dY = 0; + if( fClusterParam ){ + Int_t padSize = 0; // short pads + if (currentCluster->GetDetector() >= 36) { + padSize = 1; // medium pads + if (currentCluster->GetRow() > 63) padSize = 2; // long pads + } + dY = fClusterParam->GetWaveCorrection( padSize, + currentCluster->GetZ(), + currentCluster->GetMax(), + currentCluster->GetPad(), + riemanFitAngle.GetDYat(currentCluster->GetX()) + ); } - - //============================================================================================= - - if (useForResol && nclFound > 2 * kMinRatio * kDelta - && irow % kDeltaWriteDebugStream == 0 && fDebugLevel > 4){ - if (fDebugLevel > 5) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow); - FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta); - } // if (useForResol && nclFound > 2 * kMinRatio * kDelta) - - } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++) -} // FillResolutionHistoLocal(...) - - - -void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta) { - // - // - debug part of FillResolutionHistoLocal - - // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data - // called only for fDebugLevel > 4 - // fill resolution trees - // - - Int_t sector = cluster0->GetDetector(); - Float_t xref = cluster0->GetX(); - Int_t padSize = 0; // short pads - if (cluster0->GetDetector() >= 36) { + riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), csigmaY,csigmaZ); + riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), TMath::Abs(csigmaY*TMath::Sqrt(1+TMath::Abs(idelta))),TMath::Abs(csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta)))); + } // loop over neighbourhood for fitter filling + if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10); + if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow); + if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow + riemanFit.Update(); + riemanFitW.Update(); + Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5)); + Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5)); + Double_t paramYR[3], paramZR[3]; + Double_t paramYRW[3], paramZRW[3]; + // + paramYR[0] = riemanFit.GetYat(xref); + paramZR[0] = riemanFit.GetZat(xref); + paramYRW[0] = riemanFitW.GetYat(xref); + paramZRW[0] = riemanFitW.GetZat(xref); + // + paramYR[1] = riemanFit.GetDYat(xref); + paramZR[1] = riemanFit.GetDZat(xref); + paramYRW[1] = riemanFitW.GetDYat(xref); + paramZRW[1] = riemanFitW.GetDZat(xref); + // + Int_t reject=0; + if (chi2R>kChi2Cut) reject+=1; + if (chi2RW>kChi2Cut) reject+=2; + if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4; + if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8; + if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16; + if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32; + // + TTreeSRedirector *cstream = GetDebugStreamer(); + // get fit parameters from pol2 fit: + Double_t tracky = paramYR[0]; + Double_t trackz = paramZR[0]; + Float_t deltay = cluster0->GetY()-tracky; + Float_t deltaz = cluster0->GetZ()-trackz; + Float_t angley = paramYR[1]; + Float_t anglez = paramZR[1]; + + Int_t padSize = 0; // short pads + if (cluster0->GetDetector() >= 36) { padSize = 1; // medium pads if (cluster0->GetRow() > 63) padSize = 2; // long pads - } - - 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(); - fitZ2.ClearPoints(); - fitY2Q.ClearPoints(); - fitZ2Q.ClearPoints(); - fitY2S.ClearPoints(); - fitZ2S.ClearPoints(); - - 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; - Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) ); - Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) ); + } + Int_t ipad = TMath::Nint(cluster0->GetPad()); + if (cstream){ + Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ()); + (*cstream)<<"Resol0"<< + "run="<GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad())); + fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 ); + + + TH3F * his3 = 0; + his3 = (TH3F*)fRMSY->At(padSize); + if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) )); + his3 = (TH3F*)fRMSZ->At(padSize); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) )); + + his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) )); + his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) )); + + + his3 = (TH3F*)fResolY->At(padSize); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay ); + his3 = (TH3F*)fResolZ->At(padSize); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz ); + his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay ); + his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz ); + //============================================================================================= + // + // Fill THN histograms + // + Double_t scaleQ= TMath::Min(Double_t(cluster0->GetMax()),200.)/30.; + Bool_t isSelected = (TMath::Exp(scaleQ)>fQDownscaleRatio*gRandom->Rndm()); + if (!isSelected) continue; + Double_t xvar[9]; + xvar[1]=padSize; // pad type + xvar[2]=cluster0->GetZ(); // + xvar[3]=cluster0->GetMax(); + + xvar[0]=deltay; + double cog = cluster0->GetPad()-Int_t(cluster0->GetPad());// distance to the center of the pad + xvar[4] = cog; + xvar[5]=angley; + + if( TMath::Abs(cog-0.5)<1.e-8 ) xvar[4]=-1; // fill one-pad clusters in the underflow bin + fHisDeltaY->Fill(xvar); + + xvar[4]=cog; + xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2()); + fHisRMSY->Fill(xvar); + + xvar[0]=deltaz; + xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin()); // distance to the center of the time bin + xvar[5]=anglez; + fHisDeltaZ->Fill(xvar); + xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2()); + fHisRMSZ->Fill(xvar); + + } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++) +} // FillResolutionHistoLocal(...) -} -TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){ - // - // creates a new histogram which contains the difference between - // the histogram hfit and the function func - // - TH2D * result = (TH2D*)hfit->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() const { @@ -1268,69 +920,8 @@ void AliTPCcalibTracks::SetStyle() const { } -void AliTPCcalibTracks::Draw(Option_t* opt){ - // - // draw-function of AliTPCcalibTracks - // will draws some exemplaric pictures - // - - if (fDebugLevel > 6) Info("Draw", "Drawing an exemplaric picture."); - SetStyle(); - Double_t min = 0; - Double_t max = 0; - TCanvas *c1 = new TCanvas(); - c1->Divide(0, 3); - TVirtualPad *upperThird = c1->GetPad(1); - TVirtualPad *middleThird = c1->GetPad(2); - TVirtualPad *lowerThird = c1->GetPad(3); - upperThird->Divide(2,0); - TVirtualPad *upleft = upperThird->GetPad(1); - TVirtualPad *upright = upperThird->GetPad(2); - middleThird->Divide(2,0); - TVirtualPad *middleLeft = middleThird->GetPad(1); - TVirtualPad *middleRight = middleThird->GetPad(2); - lowerThird->Divide(2,0); - TVirtualPad *downLeft = lowerThird->GetPad(1); - TVirtualPad *downRight = lowerThird->GetPad(2); - - - upleft->cd(0); - min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20; - max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20; - fDeltaY->SetAxisRange(min, max); - fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable - c1->Update(); - - upright->cd(0); - max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20; - min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20; - fDeltaZ->SetAxisRange(min, max); - fDeltaZ->Fit("gaus","q","",min, max); - c1->Update(); - - middleLeft->cd(); - fHclus->Draw(opt); - - middleRight->cd(); - fRejectedTracksHisto->Draw(opt); - TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC"); - TText *t1 = pt->AddText("1: kEdgeThetaCutNoise"); - TText *t2 = pt->AddText("2: kMinClusters"); - TText *t3 = pt->AddText("3: kMinRatio"); - TText *t4 = pt->AddText("4: kMax1pt"); - t1 = t1; t2 = t2; t3 = t3; t4 = t4; // avoid compiler warnings - pt->SetToolTipText("Legend for failed cuts"); - pt->Draw(); - - downLeft->cd(); - fHclusterPerPadrowRaw->Draw(opt); - - downRight->cd(); - fHclusterPerPadrow->Draw(opt); -} - -void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ +void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){ // // all functions are called, that produce pictures // the histograms are written to the directory 'pathName' @@ -1338,475 +929,14 @@ void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ // 'stat' is also the number of minEntries for MakeResPlotsQTree // - if (fDebugLevel > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName); - MakeAmpPlots(stat, pathName); - MakeDeltaPlots(pathName); - FitResolutionNew(pathName); - FitRMSNew(pathName); - MakeChargeVsDriftLengthPlots(pathName); -// MakeResPlotsQ(1, 1); - MakeResPlotsQTree(stat, pathName); + if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName); + 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); - if (fDebugLevel > -1) 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); - if (fDebugLevel > -1) 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); - if (fDebugLevel > -1) cout << "creating DeltaYZ.ps..." << endl; - min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20; - max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20; - fDeltaY->SetAxisRange(min, max); - ps->NewPage(); - 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->Fit("gaus","q","",min, max); - c1->Update(); - ps->Close(); - delete ps; - delete c1; - gSystem->ChangeDirectory(".."); -} - - -void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){ - // - // creates charge vs. driftlength plots, one TProfile for each ROC - // is not correct like this, should be one TProfile for each sector and padsize - // - - SetStyle(); - gSystem->MakeDirectory(pathName); - gSystem->ChangeDirectory(pathName); - - TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable - TPostScript *ps; - ps = new TPostScript("chargeVsDriftlengthOld.ps", 112); - if (fDebugLevel > -1) cout << "creating chargeVsDriftlength.ps..." << endl; - TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone()); - TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone()); - chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC"); - chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs"); - chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC"); - chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs"); - - for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { - ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw(); - c1->Update(); - if (i > 0 && i < 36) { - chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i))); - chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36))); - } - ps->NewPage(); - } - chargeVsDriftlengthAllIROCs->Draw(); - c1->Update(); // valgrind - ps->NewPage(); - chargeVsDriftlengthAllOROCs->Draw(); - c1->Update(); - ps->Close(); - delete ps; - delete c1; - gSystem->ChangeDirectory(".."); -} -void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){ - // - // creates charge vs. driftlength plots, one TProfile for each ROC - // under development.... - // - - SetStyle(); - gSystem->MakeDirectory(pathName); - gSystem->ChangeDirectory(pathName); - - TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700)); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable -// TCanvas c1("c1", "c1", 500,(sqrt(2)*500)) - c1->Divide(0,3); - TPostScript *ps; - ps = new TPostScript("chargeVsDriftlength.ps", 111); - if (fDebugLevel > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl; - - TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone()); - TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone()); - TProfile *chargeVsDriftlengthAllLongPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone()); - chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads"); - chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads"); - chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads"); - chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads"); - chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads"); - chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads"); - - for (Int_t i = 0; i < 36; i++) { - c1->cd(1)->SetGridx(); - c1->cd(1)->SetGridy(); - ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw(); - c1->cd(2)->SetGridx(); - c1->cd(2)->SetGridy(); - ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw(); - c1->cd(3)->SetGridx(); - c1->cd(3)->SetGridy(); - ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw(); - c1->Update(); - chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)); - chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)); - chargeVsDriftlengthAllLongPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)); - ps->NewPage(); - } - c1->cd(1)->SetGridx(); - c1->cd(1)->SetGridy(); - chargeVsDriftlengthAllShortPads->Draw(); - c1->cd(2)->SetGridx(); - c1->cd(2)->SetGridy(); - chargeVsDriftlengthAllMediumPads->Draw(); - c1->cd(3)->SetGridx(); - c1->cd(3)->SetGridy(); - chargeVsDriftlengthAllLongPads->Draw(); - c1->Update(); // valgrind -// ps->NewPage(); - 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); - if (fDebugLevel > -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 * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2"); - TH3F * hisResY1 = (TH3F*)fResolY->At(1); - hisResY1->FitSlicesZ(); - TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2"); - TH3F * hisResY2 = (TH3F*)fResolY->At(2); - hisResY2->FitSlicesZ(); - TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2"); - // - ps->NewPage(); - c.cd(1); - hisResY02->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost - hisResY02->Draw("surf1"); - c.cd(2); - MakeDiff(hisResY02,fres)->Draw("surf1"); - c.Update(); - // c.SaveAs("ResolutionYPad0.eps"); - ps->NewPage(); - c.cd(1); - hisResY12->Fit(fres, "q"); - hisResY12->Draw("surf1"); - c.cd(2); - MakeDiff(hisResY12,fres)->Draw("surf1"); - c.Update(); - // c.SaveAs("ResolutionYPad1.eps"); - ps->NewPage(); - c.cd(1); - hisResY22->Fit(fres, "q"); - hisResY22->Draw("surf1"); - c.cd(2); - MakeDiff(hisResY22,fres)->Draw("surf1"); - c.Update(); -// c.SaveAs("ResolutionYPad2.eps"); - - // create histogramw for Z-resolution - TH3F * hisResZ0 = (TH3F*)fResolZ->At(0); - hisResZ0->FitSlicesZ(); - TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2"); - TH3F * hisResZ1 = (TH3F*)fResolZ->At(1); - hisResZ1->FitSlicesZ(); - TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2"); - TH3F * hisResZ2 = (TH3F*)fResolZ->At(2); - hisResZ2->FitSlicesZ(); - TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2"); - - ps->NewPage(); - c.cd(1); - hisResZ02->Fit(fres, "q"); - hisResZ02->Draw("surf1"); - c.cd(2); - MakeDiff(hisResZ02,fres)->Draw("surf1"); - c.Update(); -// c.SaveAs("ResolutionZPad0.eps"); - ps->NewPage(); - c.cd(1); - hisResZ12->Fit(fres, "q"); - hisResZ12->Draw("surf1"); - c.cd(2); - MakeDiff(hisResZ12,fres)->Draw("surf1"); - c.Update(); -// c.SaveAs("ResolutionZPad1.eps"); - ps->NewPage(); - c.cd(1); - hisResZ22->Fit(fres, "q"); - hisResZ22->Draw("surf1"); - c.cd(2); - MakeDiff(hisResZ22,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"; - hisResY02->Fit(fres, "q"); // valgrind - fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ02->Fit(fres, "q"); - fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - // - fresol<<"Pad 1.00 cm"<<1<<"\n"; - hisResY12->Fit(fres, "q"); // valgrind - fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ12->Fit(fres, "q"); - fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - // - fresol<<"Pad 1.50 cm"<<0<<"\n"; - hisResY22->Fit(fres, "q"); - fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ22->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); - if (fDebugLevel > -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 * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1"); - TH3F * hisResY1 = (TH3F*)fRMSY->At(1); - hisResY1->FitSlicesZ(); - TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1"); - TH3F * hisResY2 = (TH3F*)fRMSY->At(2); - hisResY2->FitSlicesZ(); - TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1"); - // - ps->NewPage(); - c.cd(1); - hisResY02->Fit(frms, "qn0"); - hisResY02->Draw("surf1"); - c.cd(2); - MakeDiff(hisResY02,frms)->Draw("surf1"); - c.Update(); -// c.SaveAs("RMSYPad0.eps"); - ps->NewPage(); - c.cd(1); - hisResY12->Fit(frms, "qn0"); // valgrind several blocks possibly lost - hisResY12->Draw("surf1"); - c.cd(2); - MakeDiff(hisResY12,frms)->Draw("surf1"); - c.Update(); -// c.SaveAs("RMSYPad1.eps"); - ps->NewPage(); - c.cd(1); - hisResY22->Fit(frms, "qn0"); - hisResY22->Draw("surf1"); - c.cd(2); - MakeDiff(hisResY22,frms)->Draw("surf1"); - c.Update(); -// c.SaveAs("RMSYPad2.eps"); - - // create histogramw for Z-RMS - TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0); - hisResZ0->FitSlicesZ(); - TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1"); - TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1); - hisResZ1->FitSlicesZ(); - TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1"); - TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2); - hisResZ2->FitSlicesZ(); - TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1"); - // - ps->NewPage(); - c.cd(1); - hisResZ02->Fit(frms, "qn0"); // valgrind - hisResZ02->Draw("surf1"); - c.cd(2); - MakeDiff(hisResZ02,frms)->Draw("surf1"); - c.Update(); -// c.SaveAs("RMSZPad0.eps"); - ps->NewPage(); - c.cd(1); - hisResZ12->Fit(frms, "qn0"); - hisResZ12->Draw("surf1"); - c.cd(2); - MakeDiff(hisResZ12,frms)->Draw("surf1"); - c.Update(); -// c.SaveAs("RMSZPad1.eps"); - ps->NewPage(); - c.cd(1); - hisResZ22->Fit(frms, "qn0"); // valgrind 1 block possibly lost - hisResZ22->Draw("surf1"); - c.cd(2); - MakeDiff(hisResZ22,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"; - hisResY02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost - filerms<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ02->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"; - hisResY12->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"; - hisResZ12->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"; - hisResY22->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"; - hisResZ22->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::MakeResPlotsQTree(Int_t minEntries, char* pathName){ +void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){ // // Make tree with resolution parameters // the result is written to 'resol.root' in directory 'pathname' @@ -1834,8 +964,8 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ // RMSe1: error of sigma of gaus fit in RMS-Histogram // - if (fDebugLevel > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl; - if (fDebugLevel > -1) cout << " relax, the calculation will take a while..." << endl; + if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl; + if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl; gSystem->MakeDirectory(pathName); gSystem->ChangeDirectory(pathName); @@ -1869,7 +999,7 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ } } } - if (fDebugLevel > -1) cout << "Histograms loaded, starting to proces..." << endl; + if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl; //-------------------------------------------------------------------------------------------- @@ -1889,7 +1019,7 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ // loop pad type for (Int_t iq = -1; iq < 10; iq++){ // LOOP Q - if (fDebugLevel > -1) + if (GetDebugLevel() > -1) cout << "Loop-counter, this is loop " << loopCounter << " of 66, (" << (Int_t)((loopCounter)/66.*100) << "% done), " << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush; @@ -1945,7 +1075,9 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ } qMean/=entriesQ; } - + if (!hres) continue; + if (!hrms) continue; + TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis TAxis *yAxisAngle = hres->GetYaxis(); // angle axis TAxis *zAxisDelta = hres->GetZaxis(); // delta axis @@ -1966,10 +1098,10 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ // these histograms are delta histograms for given direction, padSize, chargeBin, // angleBin and driftLengthBin // later on they will be fitted with a gausian, its sigma is the resoltuion... - sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter); + snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter); // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax()); projectionRes->SetNameTitle(name, name); - sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter); + snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter); // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax()); projectionRms->SetNameTitle(name, name); @@ -2083,7 +1215,7 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ "RMSe0="< 5) { + if (GetDebugLevel() > 5) { projectionRes->SetDirectory(fTreeResol.GetFile()); projectionRes->Write(projectionRes->GetName()); projectionRes->SetDirectory(0); @@ -2106,8 +1238,8 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ fileInfo.Write("fileInfo"); // resolFile.Close(); // fTreeResol.GetFile()->Close(); - if (fDebugLevel > -1) cout << endl; - if (fDebugLevel > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl; + if (GetDebugLevel() > -1) cout << endl; + if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl; gSystem->ChangeDirectory(".."); } @@ -2115,384 +1247,6 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ -/* -Int_t AliTPCcalibTracks::fgLoopCounter = 0; -void AliTPCcalibTracks::MakeResPlotsQTreeThread(Int_t minEntries, char* pathName){ - // - // - // - if (fDebugLevel > -1) cout << " ***** this is MakeResPlotsQTreeThread *****" << endl; - if (fDebugLevel > -1) cout << " relax, the calculation will take a while..." << endl; - if (fDebugLevel > -1) cout << " it will be done using 6 TThreads." << endl; - - gSystem->MakeDirectory(pathName); - gSystem->ChangeDirectory(pathName); - TString kFileName = "resol.root"; -// TTreeSRedirector *fTreeResol = new TTreeSRedirector(kFileName.Data()); - TTreeSRedirector fTreeResol(kFileName.Data()); - - TH3F *resArray[2][3][11]; - TH3F *rmsArray[2][3][11]; - - // load histograms from fArrayQDY and fArrayQDZ - // into resArray and rmsArray - // that is all we need here - 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); - 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); - } - } - } - if (fDebugLevel > 4) cout << "Histograms loaded, starting to proces..." << endl; - - //-------------------------------------------------------------------------------------------- - - Int_t threadCounter = 0; - TObjArray *listOfThreads = new TObjArray(); - - fgLoopCounter = 0; - for (Int_t idim = 0; idim < 2; idim++){ - // Loop y-z corrdinate - for (Int_t ipad = 0; ipad < 3; ipad++){ - // loop pad type - - // make list of variables for threads - // make threads - // list them - // execute them - - TthreadParameterStruct *structOfParameters = new TthreadParameterStruct(); - structOfParameters->logLevel = fDebugLevel; - structOfParameters->minEntries = minEntries; - structOfParameters->dim = idim; - structOfParameters->pad = ipad; - structOfParameters->resArray = &resArray; - structOfParameters->rmsArray = &rmsArray; - structOfParameters->fileName = &kFileName; - structOfParameters->fTreeResol = &fTreeResol; - TThread *thread = new TThread(Form("thread%i", threadCounter), (void(*) (void *))&MakeResPlotsQTreeThreadFunction, (void*)structOfParameters); - listOfThreads->AddAt(thread, threadCounter); - thread->Run(); -// gSystem->Sleep(500); // so that the threads do not run synchron - -// typedef TH3F test; -// test *testArray; -// TH3F* (*testArray)[2][3][11]; -// testArray = &resArray; - int i[2][3][4]; - int (*ptr)[2][3][4]; - ptr = &i; - int (*ptr2)[2][3][4] = ptr; - int j = (*ptr2)[1][1][1]; - j = j; - - threadCounter++; - } // ipad-loop - } // idim-loop - - // wait untill all threads are finished - Bool_t allFinished = kFALSE; - Int_t numberOfRunningThreads = 0; - char c[4] = {'-', '\\', '|', '/'}; - Int_t iTime = 0; - while (!allFinished) { - allFinished = kTRUE; - numberOfRunningThreads = 0; - for (Int_t i = 0; i < listOfThreads->GetEntriesFast(); i++) { - if (listOfThreads->At(i) != 0x0 && ((TThread*)listOfThreads->At(i))->GetState() == TThread::kRunningState) { - allFinished = kFALSE; - numberOfRunningThreads++; - } - } - cout << "Loop-counter, loop " << fgLoopCounter << " of 66 has just started, (" - << (Int_t)((fgLoopCounter)/66.*100) << "% done), " << "number of running TThreads: " - << numberOfRunningThreads << " \t" << c[iTime%4] << " \r" << std::flush; - iTime++; - gSystem->Sleep(500); - } - cout << endl; - - // old version: - // Real time 0:01:31, CP time 44.690 - - // version without sleep: - // Real time 0:02:18, CP time 106.280 - - // version with sleep, listOfThreads-Bug corrected: - // Real time 0:01:35, CP time 0.800 - - TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries)); - fileInfo.Write("fileInfo"); - if (fDebugLevel > -1) cout << endl; - if (fDebugLevel > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl; - gSystem->ChangeDirectory(".."); -} - - -TMutex* AliTPCcalibTracks::fgWriteMutex = new TMutex(); -TMutex* AliTPCcalibTracks::fgFitResMutex = new TMutex(); -TMutex* AliTPCcalibTracks::fgFitRmsMutex = new TMutex(); -void* AliTPCcalibTracks::MakeResPlotsQTreeThreadFunction(void* arg){ - // - // - // - - TthreadParameterStruct *structOfParameters = (TthreadParameterStruct*)arg; - Int_t fDebugLevel = structOfParameters->logLevel; - Int_t minEntries = structOfParameters->minEntries; - Int_t idim = structOfParameters->dim; - Int_t ipad = structOfParameters->pad; - TThread::Lock(); - TH3F* (*resArray)[2][3][11] = structOfParameters->resArray; - TH3F* (*rmsArray)[2][3][11] = structOfParameters->rmsArray; - TThread::UnLock(); - TString *kFileName = structOfParameters->fileName; - TTreeSRedirector *fTreeResol = structOfParameters->fTreeResol; - - if (fDebugLevel > 4) TThread::Printf("Thread started, dim = %i, pad = %i...", idim, ipad); - - TThread::Lock(); - char name[200]; - sprintf(name, "dim%ipad%i", idim, ipad); - TH1D *projectionRes = new TH1D(Form("projectionRes%s", name), "projectionRes", 50, -1, 1); - TH1D *projectionRms = new TH1D(Form("projectionRms%s", name), "projectionRms", 50, -1, 1); - char fitFuncName[200]; - sprintf(name, "fitFunctionDim%iPad%i", idim, ipad); - TF1 *fitFunction = new TF1(fitFuncName, "gaus"); - TThread::UnLock(); - - Double_t zMean, angleMean, zCenter, angleCenter; - Double_t zSigma, angleSigma; - - for (Int_t iq = -1; iq < 10; iq++){ - // LOOP Q - fgLoopCounter++; - TH3F *hres = 0; - TH3F *hrms = 0; - Double_t qMean = 0; - Float_t entriesQ = 0; - - if (fDebugLevel > 4) TThread::Printf(" start of iq-loop, dim = %i, pad = %i, iq = %i...", idim, ipad, iq); - // calculate qMean - if (iq == -1){ - // integrated spectra - 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) { - TThread::Lock(); - hres = (TH3F*)hresl->Clone(); - hrms = (TH3F*)hrmsl->Clone(); - TThread::UnLock(); - } - else{ - hres->Add(hresl); - hrms->Add(hrmsl); - } - } - qMean /= entriesQ; - qMean *= -1.; // integral mean charge - } - else { - // loop over neighboured Q-bins - // accumulate entries from neighboured Q-bins - for (Int_t iql = iq - 1; iql <= iq + 1; iql++){ - if (iql < 0) continue; - 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) { - TThread::Lock(); - hres = (TH3F*) hresl->Clone(); - hrms = (TH3F*) hrmsl->Clone(); - TThread::UnLock(); - } - else{ - hres->Add(hresl); - hrms->Add(hrmsl); - } - } - qMean/=entriesQ; - } - TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis - TAxis *yAxisAngle = hres->GetYaxis(); // angle axis - TAxis *zAxisDelta = hres->GetZaxis(); // delta axis - TAxis *zAxisRms = hrms->GetZaxis(); // rms axis - - // loop over all angle bins - for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) { - angleCenter = yAxisAngle->GetBinCenter(ibinyAngle); - // loop over all driftlength bins - for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) { - zCenter = xAxisDriftLength->GetBinCenter(ibinxDL); - zSigma = xAxisDriftLength->GetBinWidth(ibinxDL); - angleSigma = yAxisAngle->GetBinWidth(ibinyAngle); - zMean = zCenter; // changens, when more statistic is accumulated - angleMean = angleCenter; // changens, when more statistic is accumulated - - // create 2 1D-Histograms, projectionRes and projectionRms - // these histograms are delta histograms for given direction, padSize, chargeBin, - // angleBin and driftLengthBin - // later on they will be fitted with a gausian, its sigma is the resoltuion... - sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter); - projectionRes->SetNameTitle(name, name); - sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter); - projectionRms->SetNameTitle(name, name); - - projectionRes->Reset(); - projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax()); - projectionRms->Reset(); - projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax()); - projectionRes->SetDirectory(0); - projectionRms->SetDirectory(0); - - Double_t entries = 0; - Int_t nbins = 0; // counts, how many bins were accumulated - zMean = 0; - angleMean = 0; - entries = 0; - - // fill projectionRes and projectionRms for given dim, ipad and iq, - // as well as for given angleBin and driftlengthBin - // if this gives not enough statistic, include neighbourhood - // (angle and driftlength) successifely - for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin - for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction - for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction - if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time ! - Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction - Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction - if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram! - if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram! - if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram! - nbins++; // count the number of accumulated bins - // Fill resolution histo - for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) { - // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable - projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3)); - entries += hres->GetBinContent(binx2, biny2, ibin3); - zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2); - angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2); - } // ibin3 loop - // fill RMS histo - for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) { - projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3)); - } - } //dbinx2 loop - if (entries > minEntries) break; // enough statistic accumulated - } // dbiny2 loop - if (entries > minEntries) break; // enough statistic accumulated - } // dbin loop - if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree - zMean /= entries; - angleMean /= entries; - - if (entries > minEntries) { - // when enough statistic is accumulated - // fit Delta histograms with a gausian - // of the gausian is the resolution (resol), its fit error is sigma - // store also mean and RMS of the histogram - Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2; - Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2; - fitFunction->SetMaximum(xmax); - fitFunction->SetMinimum(xmin); - TThread::Lock(); - // fgFitResMutex->Lock(); - projectionRes->Fit(fitFuncName, "qN0", "", xmin, xmax); - // fgFitResMutex->UnLock(); - TThread::UnLock(); - Float_t resol = fitFunction->GetParameter(2); - Float_t sigma = fitFunction->GetParError(2); - Float_t meanR = projectionRes->GetMean(); - Float_t sigmaR = projectionRes->GetRMS(); - // fit also RMS histograms with a gausian - // store mean and sigma of the gausian in rmsMean and rmsSigma - // store also the fit errors in errorRMS and errorSigma - xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2; - xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2; - fitFunction->SetMaximum(xmax); - fitFunction->SetMinimum(xmin); - TThread::Lock(); - projectionRms->Fit(fitFuncName, "qN0", "", xmin, xmax); - TThread::UnLock(); - Float_t rmsMean = fitFunction->GetParameter(1); - Float_t rmsSigma = fitFunction->GetParameter(2); - Float_t errorRMS = fitFunction->GetParError(1); - Float_t errorSigma = fitFunction->GetParError(2); - Float_t length = 0.75; - if (ipad == 1) length = 1; - if (ipad == 2) length = 1.5; - - fgWriteMutex->Lock(); - (*fTreeResol)<<"Resol"<< - "Entries="< 5) { - projectionRes->SetDirectory(fTreeResol->GetFile()); - projectionRes->Write(projectionRes->GetName()); - projectionRes->SetDirectory(0); - projectionRms->SetDirectory(fTreeResol->GetFile()); - projectionRms->Write(projectionRms->GetName()); - projectionRes->SetDirectory(0); - } - fgWriteMutex->UnLock(); - } // if (projectionRes->GetSum() > minEntries) - } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) - } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) - - } // iq-loop - delete projectionRes; - delete projectionRms; - if (fDebugLevel > 4) TThread::Printf("Ende, dim = %i, pad = %i", idim, ipad); - return 0; -} - -*/ - - Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { // // function to merge several AliTPCcalibTracks objects after PROOF calculation @@ -2500,12 +1254,12 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!! // - if (fDebugLevel > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl; + if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl; if (!collectionList) return 0; if (collectionList->IsEmpty()) return -1; - if (fDebugLevel > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1 - if (fDebugLevel > 5) cout << " the list in the merge-function looks as follows: " << endl; + if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1 + if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl; collectionList->Print(); // create a list for each data member @@ -2513,17 +1267,12 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { TList* deltaZList = new TList; TList* arrayAmpRowList = new TList; TList* rejectedTracksList = new TList; - TList* hclusList = new TList; TList* clusterCutHistoList = new TList; TList* arrayAmpList = new TList; TList* arrayQDYList = new TList; TList* arrayQDZList = new TList; TList* arrayQRMSYList = new TList; TList* arrayQRMSZList = new TList; - TList* arrayChargeVsDriftlengthList = new TList; - TList* calPadRegionChargeVsDriftlengthList = new TList; - TList* hclusterPerPadrowList = new TList; - TList* hclusterPerPadrowRawList = new TList; TList* resolYList = new TList; TList* resolZList = new TList; TList* rMSYList = new TList; @@ -2535,16 +1284,11 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { TIterator *listIterator = collectionList->MakeIterator(); AliTPCcalibTracks *calibTracks = 0; - if (fDebugLevel > 1) cout << "start to iterate, filling lists" << endl; + if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl; Int_t counter = 0; - while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){ + while ( (calibTracks = dynamic_cast (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()); @@ -2553,75 +1297,36 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { resolZList->Add(calibTracks->GetfResolZ()); rMSYList->Add(calibTracks->GetfRMSY()); rMSZList->Add(calibTracks->GetfRMSZ()); - arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength()); - calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength()); - hclusList->Add(calibTracks->GetfHclus()); rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto()); clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto()); - hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow()); - hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw()); - fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad()); - fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw()); + // + if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad()) + fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad()); + // fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw()); counter++; - if (fDebugLevel > 5) cout << "filling lists, object " << counter << " added." << endl; + if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl; + AddHistos(calibTracks); } // merge data members - if (fDebugLevel > 0) cout << "histogram's merge-functins are called... " << endl; - fDeltaY->Merge(deltaYList); - fDeltaZ->Merge(deltaZList); - fHclus->Merge(hclusList); + if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl; fClusterCutHisto->Merge(clusterCutHistoList); fRejectedTracksHisto->Merge(rejectedTracksList); - fHclusterPerPadrow->Merge(hclusterPerPadrowList); - fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList); TObjArray* objarray = 0; TH1* hist = 0; TList* histList = 0; TIterator *objListIterator = 0; - if (fDebugLevel > 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; - } - - if (fDebugLevel > 0) 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; - } - - if (fDebugLevel > 0) cout << "merging fArrayQDY..." << endl; + + if (GetDebugLevel() > 0) 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); } @@ -2630,14 +1335,13 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete objListIterator; } - if (fDebugLevel > 0) cout << "merging fArrayQDZ..." << endl; + if (GetDebugLevel() > 0) 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); } @@ -2646,14 +1350,14 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete objListIterator; } - if (fDebugLevel > 0) cout << "merging fArrayQRMSY..." << endl; + if (GetDebugLevel() > 0) 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; + // if (!objarray) continue; // removed for coverity -> JMT hist = (TH3F*)(objarray->At(i)); histList->Add(hist); } @@ -2662,14 +1366,13 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete objListIterator; } - if (fDebugLevel > 0) cout << "merging fArrayQRMSZ..." << endl; + if (GetDebugLevel() > 0) 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); } @@ -2678,65 +1381,18 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete objListIterator; } - if (fDebugLevel > 0) cout << "merging fArrayChargeVsDriftlength..." << endl; - // merge fArrayChargeVsDriftlength - for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300 - objListIterator = arrayChargeVsDriftlengthList->MakeIterator(); - histList = new TList; - while (( objarray = (TObjArray*)objListIterator->Next() )) { - // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile - if (!objarray) continue; - hist = (TProfile*)(objarray->At(i)); - histList->Add(hist); - } - ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList); - delete histList; - delete objListIterator; - } - - if (fDebugLevel > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl; - // merge fcalPadRegionChargeVsDriftlength - AliTPCCalPadRegion *cpr = 0x0; - - /* - TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator(); - while (hist = (TProfile*)regionIterator->Next()) { - // loop over all calPadRegion's in destination calibTracks object - objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator(); - while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) { - - - hist->Merge(...); - } - */ - for (UInt_t isec = 0; isec < 36; isec++) { - for (UInt_t padSize = 0; padSize < 3; padSize++){ - objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator(); - histList = new TList; - while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) { - // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object - if (!cpr) continue; - hist = (TProfile*)cpr->GetObject(isec, padSize); - histList->Add(hist); - } - ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList); - delete histList; - delete objListIterator; - } - } - if (fDebugLevel > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl; + if (GetDebugLevel() > 0) 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); } @@ -2751,7 +1407,6 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { 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); } @@ -2766,7 +1421,6 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { 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); } @@ -2781,7 +1435,6 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { 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); } @@ -2802,70 +1455,659 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete resolZList; delete rMSYList; delete rMSZList; -// delete nRowsList; -// delete nSectList; -// delete fileNoList; delete listIterator; - if (fDebugLevel > 0) cout << "merging done!" << endl; + if (GetDebugLevel() > 0) 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 - .L libTPCcalib.so - 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++) { - if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl; - list->Add(new AliTPCcalibTracks(ct)); - } - - // 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; + +void AliTPCcalibTracks::MakeHistos(){ + // + ////make THnSparse + // + //THnSparse *fHisDeltaY; // THnSparse - delta Y + //THnSparse *fHisDeltaZ; // THnSparse - delta Z + //THnSparse *fHisRMSY; // THnSparse - rms Y + //THnSparse *fHisRMSZ; // THnSparse - rms Z + //THnSparse *fHisQmax; // THnSparse - qmax + //THnSparse *fHisQtot; // THnSparse - qtot + // cluster performance bins + // 0 - variable of interest + // 1 - pad type - 0- short 1-medium 2-long pads + // 2 - drift length - drift length -0-1 + // 3 - Qmax - Qmax - 2- 400 + // 4 - cog - COG position - 0-1 + // 5 - tan(phi) - local y angle + // 6 - tan(theta) - local z angle + // 7 - sector - sector number + Double_t xminTrack[8], xmaxTrack[8]; + Int_t binsTrack[8]; + TString axisName[8]; - - AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018); - AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5); - 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; + // + binsTrack[0]=200; + axisName[0] ="var"; + + binsTrack[1] =3; + xminTrack[1] =0; xmaxTrack[1]=3; + axisName[1] ="pad type"; + // + binsTrack[2] =20; + xminTrack[2] =-250; xmaxTrack[2]=250; + axisName[2] ="z"; + // + binsTrack[3] =10; + xminTrack[3] =1; xmaxTrack[3]=400; + axisName[3] ="Qmax"; + // + binsTrack[4] =20; + xminTrack[4] =0; xmaxTrack[4]=1; + axisName[4] ="cog"; + // + binsTrack[5] =15; + xminTrack[5] =-1.5; xmaxTrack[5]=1.5; + axisName[5] ="tan(angle)"; + // + // + xminTrack[0] =-1.5; xmaxTrack[0]=1.5; + fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =-1.5; xmaxTrack[0]=1.5; + fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =0.; xmaxTrack[0]=0.5; + fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =0.; xmaxTrack[0]=0.5; + fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =0.; xmaxTrack[0]=100; + fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack); + + xminTrack[0] =0.; xmaxTrack[0]=250; + fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack); + + + for (Int_t ivar=0;ivar<6;ivar++){ + fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data()); + fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data()); + fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data()); + } + + + BinLogX(fHisDeltaY,3); + BinLogX(fHisDeltaZ,3); + BinLogX(fHisRMSY,3); + BinLogX(fHisRMSZ,3); + BinLogX(fHisQmax,3); + BinLogX(fHisQtot,3); + +} +void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){ + // + // Add histograms + // + if (!calib->fHisDeltaY) return; + if (calib->fHisDeltaY->GetEntries()> fgkMergeEntriesCut) return; + if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY); + if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ); + if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY); + if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ); } +void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){ + // + // Dump summary info + // + // 0.OBJ: TAxis var + // 1.OBJ: TAxis pad type + // 2.OBJ: TAxis z + // 3.OBJ: TAxis Qmax + // 4.OBJ: TAxis cog + // 5.OBJ: TAxis tan(angle) + // + if (ptype>3) return; + Int_t idim[6]={0,1,2,3,4,5}; + TString hname[4]={"dy","dz","rmsy","rmsz"}; + // + Int_t nbins5=hisInput->GetAxis(5)->GetNbins(); + Int_t first5=hisInput->GetAxis(5)->GetFirst(); + Int_t last5 =hisInput->GetAxis(5)->GetLast(); + // + for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){ // axis 5 - local angle + Bool_t bin5All=kTRUE; + if (ibin5>=first5){ + hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5)); + if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5)); + bin5All=kFALSE; + } + Double_t x5= hisInput->GetAxis(5)->GetBinCenter(ibin5); + THnSparse * his5= hisInput->Projection(5,idim); //projected histogram according selection 5 + Int_t nbins4=his5->GetAxis(4)->GetNbins(); + Int_t first4=his5->GetAxis(4)->GetFirst(); + Int_t last4 =his5->GetAxis(4)->GetLast(); + + for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){ // axis 4 - cog + Bool_t bin4All=kTRUE; + if (ibin4>=first4){ + his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4)); + if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4)); + bin4All=kFALSE; + } + Double_t x4= his5->GetAxis(4)->GetBinCenter(ibin4); + THnSparse * his4= his5->Projection(4,idim); //projected histogram according selection 4 + // + Int_t nbins3=his4->GetAxis(3)->GetNbins(); + Int_t first3=his4->GetAxis(3)->GetFirst(); + Int_t last3 =his4->GetAxis(3)->GetLast(); + // + for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - Qmax + Bool_t bin3All=kTRUE; + if (ibin3>=first3){ + his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3)); + bin3All=kFALSE; + } + Double_t x3= his4->GetAxis(3)->GetBinCenter(ibin3); + THnSparse * his3= his4->Projection(3,idim); //projected histogram according selection 3 + // + Int_t nbins2 = his3->GetAxis(2)->GetNbins(); + Int_t first2 = his3->GetAxis(2)->GetFirst(); + Int_t last2 = his3->GetAxis(2)->GetLast(); + // + for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){ // axis 2 - z + Bool_t bin2All=kTRUE; + Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2); + if (ibin2>=first2){ + his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2)); + if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2)); + bin2All=kFALSE; + } + THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2 + // + Int_t nbins1 = his2->GetAxis(1)->GetNbins(); + Int_t first1 = his2->GetAxis(1)->GetFirst(); + Int_t last1 = his2->GetAxis(1)->GetLast(); + for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){ //axis 1 - pad type + Bool_t bin1All=kTRUE; + if (ibin1>=first1){ + his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1)); + bin1All=kFALSE; + } + Double_t x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5); + TH1 * hisDelta = his2->Projection(0); + Double_t entries = hisDelta->GetEntries(); + Double_t mean=0, rms=0; + if (entries>10){ + mean = hisDelta->GetMean(); + rms = hisDelta->GetRMS(); + hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms); + mean = hisDelta->GetMean(); + rms = hisDelta->GetRMS(); + } + // + (*pcstream)<GetNdimensions(); + Long64_t nFilledBins = H->GetNbins(); + Long64_t nStatBins = 0; + + Int_t *nBins = new Int_t [nDim]; + Double_t *xMin = new Double_t [nDim]; + Double_t *xMax = new Double_t [nDim]; + Int_t *idx = new Int_t [nDim]; + + TString nameMean = H->GetName(); + TString nameSigma = H->GetName(); + TString nameEntr = H->GetName(); + nameMean+="_Mean"; + nameSigma+="_Sigma"; + nameEntr+="_Entr"; + + ok = ok && ( nBins && xMin && xMax && idx ); + + if( ok ){ + for( int i=0; iGetAxis(i)->GetXmin(); + xMax[i] = H->GetAxis(i)->GetXmax(); + nBins[i] = H->GetAxis(i)->GetNbins(); + } + + Mean = new THnSparseF( nameMean.Data(), nameMean.Data(), nDim-1, nBins+1, xMin+1, xMax+1 ); + Sigma = new THnSparseF( nameSigma.Data(), nameSigma.Data(), nDim-1, nBins+1, xMin+1, xMax+1 ); + Entr = new THnSparseF( nameEntr.Data(), nameEntr.Data(), nDim-1, nBins+1, xMin+1, xMax+1 ); + } + + ok = ok && ( Mean && Sigma && Entr ); + + if( ok ){ + for( int i=0; iGetAxis(i+1); + Mean->GetAxis(i)->SetName(ax->GetName()); + Mean->GetAxis(i)->SetTitle(ax->GetTitle()); + Sigma->GetAxis(i)->SetName(ax->GetName()); + Sigma->GetAxis(i)->SetTitle(ax->GetTitle()); + Entr->GetAxis(i)->SetName(ax->GetName()); + Entr->GetAxis(i)->SetTitle(ax->GetTitle()); + if( ax->GetXbins()->fN>0 ){ + Mean->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray()); + Sigma->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray()); + Entr->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray()); + } + } + + // Allocate bins + + for( Long64_t binS=0; binSGetBinContent(binS,idx); + Mean->GetBin(idx+1,kTRUE); + Sigma->GetBin(idx+1,kTRUE); + Entr->GetBin(idx+1,kTRUE); + } + + nStatBins = Mean->GetNbins(); + + } + + + Long64_t *hMap = new Long64_t[nFilledBins]; + Double_t *hVal = new Double_t[nFilledBins]; + Double_t *hEntr = new Double_t[nFilledBins]; + Double_t *meanD = new Double_t[nStatBins]; + Double_t *sigmaD = new Double_t[nStatBins]; + + ok = ok && ( hMap && hVal && hEntr && meanD && sigmaD ); + + // Map bins to Stat bins + + if( ok ){ + for( Long64_t binS=0; binSGetBinContent(binS,idx); + Long64_t bin = Mean->GetBin(idx+1,kFALSE); + ok = ok && ( bin>=0 && binGetBin(idx+1,kFALSE) && bin== Entr->GetBin(idx+1,kFALSE) ); + if( !ok ) break; + if( val < 0 ){ + cout << "AliTPCcalibTracks: GetSparseStat : Unexpected content of the input THn"<nBins[0] ) bin = -1; + hMap[binS] = bin; + hVal[binS] = idx[0]; + hEntr[binS] = val; + } + } + + // do 2 iteration with cut at 4 Sigma + + for( int iter=0; ok && iter<2; iter++ ){ + + // clean up entries + + for( Long64_t bin=0; bin < nStatBins; bin++){ + Entr->SetBinContent(bin, 0); + meanD[bin]=0; + sigmaD[bin]=0; + } + + // get Entries and Mean + + for( Long64_t binS=0; binSGetBinContent(bin); + double d = val - Mean->GetBinContent(bin); + if( d*d>s2*16 ) continue; + } + meanD[bin]+= entr*val; + Entr->AddBinContent(bin,entr); + } + + for( Long64_t bin=0; binGetBinContent(bin); + if( sum>=1 ){ + Mean->SetBinContent(bin, val/sum ); + } + else Mean->SetBinContent(bin, 0); + Entr->SetBinContent(bin, 0); + } + + // get RMS + + for( Long64_t binS=0; binSGetBinContent(bin); + d2 *= d2; + if( iter!=0 ){ // cut + double s2 = Sigma->GetBinContent(bin); + if( d2>s2*16 ) continue; + } + sigmaD[bin] += entr*d2; + Entr->AddBinContent(bin,entr); + } + + for( Long64_t bin=0; binGetBinContent(bin); + if( sum>=1 && val>=0 ){ + Sigma->SetBinContent(bin, val/sum ); + } + else Sigma->SetBinContent(bin, 0); + } + } + + // scale the Mean and the Sigma + + if( ok ){ + double v0 = H->GetAxis(0)->GetBinCenter(1); + double vscale = H->GetAxis(0)->GetBinWidth(1); + + for( Long64_t bin=0; binGetBinContent(bin); + double s2 = Sigma->GetBinContent(bin); + double sum = Entr->GetBinContent(bin); + if( sum>0 && s2>=0 ){ + Mean->SetBinContent(bin, v0 + (m-1)*vscale ); + Sigma->SetBinContent(bin, sqrt(s2)*vscale ); + } + else{ + Mean->SetBinContent(bin, 0); + Sigma->SetBinContent(bin, 0); + Entr->SetBinContent(bin, 0); + } + } + } + + delete[] nBins; + delete[] xMin; + delete[] xMax; + delete[] idx; + delete[] hMap; + delete[] hVal; + delete[] hEntr; + delete[] meanD; + delete[] sigmaD; + + if( !ok ){ + cout << "AliTPCcalibTracks: GetSparseStat : No memory or internal Error "<GetNdimensions(); + if( nDim<6 ){ + cout << "AliTPCcalibTracks: CreateWaveCorrection: Unknown input"<GetAxis(i)->GetXmin(); + xMax[i] = DeltaY->GetAxis(i)->GetXmax(); + nBins[i] = DeltaY->GetAxis(i)->GetNbins(); + nBinsNew[i] = nBins[i]; + } + + // Merge cog axis + if( MirrorPad ){ + Int_t centralBin = DeltaY->GetAxis(4)->FindFixBin(0.5); + xMin[4] = DeltaY->GetAxis(4)->GetBinLowEdge(centralBin); + nBinsNew[4] = nBinsNew[4]-centralBin+1; + } + + // Merge Z axis + if( MirrorZ ){ + Int_t centralBin = DeltaY->GetAxis(2)->FindFixBin(0.0); + xMin[2] = DeltaY->GetAxis(2)->GetBinLowEdge(centralBin)-0.0; + nBinsNew[2] = nBinsNew[2]-centralBin+1; + } + + // Merge Angle axis + if( MirrorAngle ){ + Int_t centralBin = DeltaY->GetAxis(5)->FindFixBin(0.0); + xMin[5] = DeltaY->GetAxis(5)->GetBinLowEdge(centralBin)-0.0; + nBinsNew[5] = nBinsNew[5]-centralBin+1; + } + + // Merge Sparse + + mergedDeltaY = new THnSparseD("mergedDeltaY", "mergedDeltaY",nDim, nBinsNew, xMin, xMax ); + + if( !mergedDeltaY ){ + cout << "AliTPCcalibTracks: CreateWaveCorrection: Can not copy a Sparse"<GetAxis(i); + mergedDeltaY->GetAxis(i)->SetName(ax->GetName()); + mergedDeltaY->GetAxis(i)->SetTitle(ax->GetTitle()); + if( ax->GetXbins()->fN>0 ){ + mergedDeltaY->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray()); + } + } + + for( Long64_t binS=0; binSGetNbins(); binS++){ + Double_t stat = DeltaY->GetBinContent(binS,idx); + if( stat < 1 ) continue; + bool swap0=0; + + if( MirrorPad && idx[4]>0){ // underflow reserved for contains one-pad clusters, don't mirror + Double_t v = DeltaY->GetAxis(4)->GetBinCenter(idx[4]); + if( v < 0.5 ) swap0 = !swap0; + idx[4] = mergedDeltaY->GetAxis(4)->FindFixBin( 0.5 + TMath::Abs(0.5 - v) ); + } + + if( MirrorZ ){ + Double_t v = DeltaY->GetAxis(2)->GetBinCenter(idx[2]); + if( v < 0.0 ) swap0 = !swap0; + if( idx[2]<=0 ) idx[2] = nBinsNew[2]+1; + else idx[2] = mergedDeltaY->GetAxis(2)->FindFixBin( TMath::Abs(v) ); + } + + if( MirrorAngle ){ + Double_t v = DeltaY->GetAxis(5)->GetBinCenter(idx[5]); + if( idx[5]<=0 ) idx[5] = nBinsNew[5]+1; + else idx[5] = mergedDeltaY->GetAxis(5)->FindFixBin( TMath::Abs(v) ); + } + + if( swap0 ){ + if( idx[0]<=0 ) idx[0] = nBinsNew[0]+1; + else if( idx[0] >= nBins[0]+1 ) idx[0] = 0; + else { + Double_t v = DeltaY->GetAxis(0)->GetBinCenter(idx[0]); + idx[0] = mergedDeltaY->GetAxis(0)->FindFixBin(-v); + } + } + + Long64_t bin = mergedDeltaY->GetBin(idx,kTRUE); + if( bin<0 ){ + cout << "AliTPCcalibTracks: CreateWaveCorrection : wrong bining"<AddBinContent(bin,stat); + } + + ret = GetTHnStat(mergedDeltaY, MeanY, SigmaY, EntrY ); + } + + if( ret==0 ){ + + MeanY->SetName("TPCWaveCorrectionMap"); + MeanY->SetTitle("TPCWaveCorrectionMap"); + SigmaY->SetName("TPCResolutionMap"); + SigmaY->SetTitle("TPCResolutionMap"); + EntrY->SetName("TPCWaveCorrectionEntr"); + EntrY->SetTitle("TPCWaveCorrectionEntr"); + + for( Long64_t bin=0; binGetNbins(); bin++){ + Double_t stat = EntrY->GetBinContent(bin,idx); + + // Normalize : Set no correction for one pad clusters + + if( idx[3]<=0 ) MeanY->SetBinContent(bin,0); + + // Suppress bins with low statistic + + if( statSetBinContent(bin,0); + MeanY->SetBinContent(bin,0); + SigmaY->SetBinContent(bin,-1); + } + } + + } + + delete[] nBins; + delete[] nBinsNew; + delete[] xMin; + delete[] xMax; + delete[] idx; + delete mergedDeltaY; + + if( ret!=0 ){ + delete MeanY; + delete SigmaY; + delete EntrY; + MeanY = 0; + SigmaY = 0; + EntrY = 0; + } + + return ret; +} + +int AliTPCcalibTracks::UpdateClusterParam( AliTPCClusterParam *cParam, Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat ) +{ + if( !cParam || !fHisDeltaY) return -1; + + THnBase *meanY = 0; + THnBase *sigmaY = 0; + THnBase *entrY = 0; + int ret = CreateWaveCorrection(fHisDeltaY, meanY, sigmaY, entrY, MirrorZ, MirrorAngle, MirrorPad, MinStat ); + if( ret<0 ) return ret; + cParam->SetWaveCorrectionMap(meanY); + cParam->SetResolutionYMap(sigmaY); + delete meanY; + delete sigmaY; + delete entrY; + return 0; +}