X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCcalibTracks.cxx;h=39d8f4919eec5828bbbbecc6dd1b025c0ea859a0;hb=8c6128b71eea83b6ee9f8bfe1de1c4f37915c233;hp=700bfaef615b99eb545dabdeeeb724b513947ea3;hpb=46e897933a84394e980cc941064170a9d61dbfac;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCcalibTracks.cxx b/TPC/AliTPCcalibTracks.cxx index 700bfaef615..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,23 @@ 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" + +*/ // /////////////////////////////////////////////////////////////////////////////// @@ -98,8 +152,10 @@ using namespace std; #include "TCut.h" #include "THnSparse.h" #include "AliRieman.h" +#include "TRandom.h" +Double_t AliTPCcalibTracks::fgkMergeEntriesCut=10000000.; //10**7 clusters ClassImp(AliTPCcalibTracks) @@ -114,6 +170,8 @@ AliTPCcalibTracks::AliTPCcalibTracks(): 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), @@ -146,6 +204,8 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): 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), @@ -230,6 +290,8 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al 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), @@ -413,6 +475,10 @@ void AliTPCcalibTracks::Process(AliTPCseed *track){ // FillResolutionHistoLocal(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) { @@ -523,6 +589,18 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ // and to avoid redundant data // + /* {//SG + static long n=0; + if( n==0 ){ + if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){ + cout<<"Map found "<GetWaveCorrectionMap() ){ + for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ + // 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( riemanFitAngle.GetN()>3 ) riemanFitAngle.Update(); + } + + // do fit + AliRieman riemanFit(2*kDelta); AliRieman riemanFitW(2*kDelta); for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ @@ -640,8 +739,23 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ if (idelta > 0){ ncl1++; } - riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ); - riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY*TMath::Sqrt(1+TMath::Abs(idelta)),csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta))); + //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()) + ); + } + 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); @@ -733,14 +847,14 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ TH3F * his3 = 0; his3 = (TH3F*)fRMSY->At(padSize); - if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) ); + if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) )); his3 = (TH3F*)fRMSZ->At(padSize); - if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) ); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(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(cluster0->GetSigmaY2()) ); + 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(cluster0->GetSigmaZ2()) ); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) )); his3 = (TH3F*)fResolY->At(padSize); @@ -755,20 +869,28 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ // // 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; - xvar[2]=cluster0->GetZ(); + xvar[1]=padSize; // pad type + xvar[2]=cluster0->GetZ(); // xvar[3]=cluster0->GetMax(); xvar[0]=deltay; - xvar[4]=cluster0->GetPad()-Int_t(cluster0->GetPad()); - xvar[5]=angley; - fHisDeltaY->Fill(xvar); + 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()); + 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()); @@ -976,10 +1098,10 @@ void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const 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); @@ -1220,7 +1342,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); } @@ -1236,7 +1357,7 @@ 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; + // if (!objarray) continue; // removed for coverity -> JMT hist = (TH3F*)(objarray->At(i)); histList->Add(hist); } @@ -1252,7 +1373,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); } @@ -1273,7 +1393,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); } @@ -1288,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); } @@ -1303,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); } @@ -1318,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); } @@ -1441,6 +1557,8 @@ 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); @@ -1541,17 +1659,17 @@ void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *p // (*pcstream)<GetAxis(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; +}