X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCcalibTracks.cxx;h=d89e22632bf622dfcd00071712f4a824bf437b96;hb=7e124979ac4d7074184f476e15f24634f2f4183f;hp=ea5582d78327158f3860a87fb14ec2771a4f9e9f;hpb=c32da87924a459f46841c6e2f09acf834b1db1a2;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCcalibTracks.cxx b/TPC/AliTPCcalibTracks.cxx index ea5582d7832..d89e22632bf 100644 --- a/TPC/AliTPCcalibTracks.cxx +++ b/TPC/AliTPCcalibTracks.cxx @@ -34,6 +34,26 @@ Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam) */ +/* + +How to retrive it from file (created using calibration task): + +gSystem->Load("libANALYSIS"); +gSystem->Load("libTPCcalib"); +TFile fcalib("CalibObjects.root"); +TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); +AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks"); + + +//USAGE of debug stream example + gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); + gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+") + AliXRDPROOFtoolkit tool; + TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200); + chainres->Lookup(); +*/ + + // /////////////////////////////////////////////////////////////////////////////// @@ -48,6 +68,8 @@ using namespace std; #include #include #include +#include + // //#include #include @@ -66,6 +88,7 @@ using namespace std; #include #include #include +#include // // AliROOT includes @@ -91,11 +114,9 @@ using namespace std; #include "TText.h" #include "TPaveText.h" #include "TSystem.h" +#include "TStatToolkit.h" +#include "TCut.h" -// Thread-stuff -//#include "TThread.h" -//#include "TMutex.h" -//#include "TLockFile.h" ClassImp(AliTPCcalibTracks) @@ -104,7 +125,6 @@ ClassImp(AliTPCcalibTracks) AliTPCcalibTracks::AliTPCcalibTracks(): AliTPCcalibBase(), fClusterParam(0), - fDebugStream(0), fROC(0), fArrayAmpRow(0), fArrayAmp(0), @@ -127,27 +147,20 @@ AliTPCcalibTracks::AliTPCcalibTracks(): 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; + SetDebugLevel(1); + if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl; } AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): - AliTPCcalibBase(), + AliTPCcalibBase(calibTracks), fClusterParam(0), - fDebugStream(0), fROC(0), fArrayAmpRow(0), fArrayAmp(0), @@ -170,19 +183,12 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): 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); @@ -240,7 +246,6 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(), calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise()); - fDebugLevel = calibTracks.GetLogLevel(); SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle()); TH1::AddDirectory(dirStatus); // set status back to original status // cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED @@ -262,7 +267,6 @@ AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibT AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) : AliTPCcalibBase(), fClusterParam(0), - fDebugStream(0), fROC(0), fArrayAmpRow(0), fArrayAmp(0), @@ -285,14 +289,7 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al 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 @@ -300,14 +297,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. // this->SetName(name); this->SetTitle(title); - if (fDebugLevel > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl; + if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl; G__SetCatchException(0); fClusterParam = clusterParam; @@ -318,8 +316,7 @@ 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); TH1::AddDirectory(kFALSE); @@ -327,7 +324,7 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al 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); + fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -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); fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad"); @@ -353,22 +350,22 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al // amplitude sprintf(chname,"Amp_Sector%d",i); - his1 = new TH1F(chname,chname,250,0,500); // valgrind + his1 = new TH1F(chname,chname,100,0,32); // 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 = new TH1F(chname,chname,100,0,32); // 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 = new TProfile(chname, chname, 25, 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 = new TProfile(chname, chname, 25, 0, 250); prof1->SetYTitle("Charge"); prof1->SetXTitle("Driftlength"); fArrayChargeVsDriftlength->AddAt(prof1,i+36); @@ -424,18 +421,18 @@ 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]; + sprintf(hname,"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); + sprintf(hname,"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); + sprintf(hname,"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); + sprintf(hname,"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); } } @@ -454,14 +451,8 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al } } - 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 } @@ -471,7 +462,7 @@ 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++){ @@ -509,13 +500,7 @@ AliTPCcalibTracks::~AliTPCcalibTracks() { delete fArrayChargeVsDriftlength->At(i); } - delete fFitterLinY1; - delete fFitterLinZ1; - delete fFitterLinY2; - delete fFitterLinZ2; - delete fFitterParY; - delete fFitterParZ; - + delete fArrayQDY; delete fArrayQDZ; delete fArrayQRMSY; @@ -529,22 +514,13 @@ AliTPCcalibTracks::~AliTPCcalibTracks() { delete fHclusterPerPadrowRaw; if (fCalPadClusterPerPad) delete fCalPadClusterPerPad; if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw; - fcalPadRegionChargeVsDriftlength->Delete(); - delete fcalPadRegionChargeVsDriftlength; - if (fDebugLevel > 4) delete fDebugStream; + if(fcalPadRegionChargeVsDriftlength) { + fcalPadRegionChargeVsDriftlength->Delete(); + delete fcalPadRegionChargeVsDriftlength; + } } -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){ // @@ -553,7 +529,7 @@ void AliTPCcalibTracks::Process(AliTPCseed *track){ // FillResolutionHistoLocal(track) // AlignUpDown(track, esd) // - if (fDebugLevel > 5) Info("Process","Starting to process the track..."); + if (GetDebugLevel() > 5) Info("Process","Starting to process the track..."); Int_t accpetStatus = AcceptTrack(track); if (accpetStatus == 0) { FillResolutionHistoLocal(track); @@ -634,7 +610,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; } @@ -666,20 +642,29 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ // and to avoid redundant data // - if (fDebugLevel > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****"); + static TLinearFitter fFitterLinY1(2,"pol1"); // + static TLinearFitter fFitterLinZ1(2,"pol1"); // + static TLinearFitter fFitterLinY2(2,"pol1"); // + static TLinearFitter fFitterLinZ2(2,"pol1"); // + static TLinearFitter fFitterParY(3,"pol2"); // + static TLinearFitter fFitterParZ(3,"pol2"); // + + fFitterLinY1.StoreData(kFALSE); + fFitterLinZ1.StoreData(kFALSE); + fFitterLinY2.StoreData(kFALSE); + fFitterLinZ2.StoreData(kFALSE); + fFitterParY.StoreData(kFALSE); + fFitterParZ.StoreData(kFALSE); + + + if (GetDebugLevel() > 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 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); @@ -716,27 +701,27 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ if (sector != sectorG){ // track leaves sector before it crossed enough rows to fit / initialization nClusters = 0; - fFitterParY->ClearPoints(); - fFitterParZ->ClearPoints(); + fFitterParY.ClearPoints(); + fFitterParZ.ClearPoints(); sectorG = sector; } else { nClusters++; Double_t x = cluster0->GetX(); - fFitterParY->AddPoint(&x, cluster0->GetY(), 1); - fFitterParZ->AddPoint(&x, cluster0->GetZ(), 1); + 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(); + fFitterParY.Eval(); + fFitterParZ.Eval(); nTrackletsAll++; - csigmaY += fFitterParY->GetChisquare() / (nClusters - 3.); - csigmaZ += fFitterParZ->GetChisquare() / (nClusters - 3.); + csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.); + csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.); nClusters = -1; - fFitterParY->ClearPoints(); - fFitterParZ->ClearPoints(); + fFitterParY.ClearPoints(); + fFitterParZ.ClearPoints(); } } } // for (Int_t irow = 0; irow < 159; irow++) @@ -758,12 +743,12 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ Float_t xref = cluster0->GetX(); // Make Fit - fFitterParY->ClearPoints(); - fFitterParZ->ClearPoints(); - fFitterLinY1->ClearPoints(); - fFitterLinZ1->ClearPoints(); - fFitterLinY2->ClearPoints(); - fFitterLinZ2->ClearPoints(); + 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 @@ -786,89 +771,102 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ nclFound++; if (idelta < 0){ ncl0++; - fFitterLinY1->AddPoint(&x, currentCluster->GetY(), csigmaY); - fFitterLinZ1->AddPoint(&x, currentCluster->GetZ(), csigmaZ); + 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); + fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY); + fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ); } - fFitterParY->AddPoint(&x, currentCluster->GetY(), csigmaY); - fFitterParZ->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 - + 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 + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + (*cstream)<<"Cut9"<< + "chi2="< 4){ - fFitterLinY1->Eval(); - fFitterLinZ1->Eval(); + fFitterLinY1.Eval(); + fFitterLinZ1.Eval(); } if (ncl1 > 4){ - fFitterLinY2->Eval(); - fFitterLinZ2->Eval(); + fFitterLinY2.Eval(); + fFitterLinZ2.Eval(); } 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); + 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; + Double_t cchi2 = 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); + cchi2 += 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); + cchi2 += 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 + // REMOVE KINK - TO be fixed - proper chi2 calculation for curved track to be implemented + //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 (cstream){ + (*cstream)<<"Cut8"<< + "chi2="<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); + 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]; @@ -891,16 +889,17 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector); if ( irow >= kFirstLargePad) max /= kLargePadSize; - profAmpRow->Fill( (Double_t)cluster0->GetRow(), max ); + Double_t smax = TMath::Sqrt(max); + profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax ); TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector); - hisAmp->Fill(max); + hisAmp->Fill(smax); // remove the following two lines one day: TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector); - profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max ); + profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax ); TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize)); - profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max ); + profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax ); fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow Int_t ipad = TMath::Nint(cluster0->GetPad()); @@ -922,8 +921,29 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ // Fill resolution histograms Bool_t useForResol = kTRUE; - if (fFitterParY->GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE; - + if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE; + + if (cstream){ + Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ()); + Float_t sy = cluster0->GetSigmaY2(); + Float_t sz = cluster0->GetSigmaZ2(); + (*cstream)<<"Resol0"<< + "run="<Fill(deltay); fDeltaZ->Fill(deltaz); @@ -940,8 +960,8 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ //============================================================================================= if (useForResol && nclFound > 2 * kMinRatio * kDelta - && irow % kDeltaWriteDebugStream == 0 && fDebugLevel > 4){ - if (fDebugLevel > 5) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow); + && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){ + if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow); FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta); } // if (useForResol && nclFound > 2 * kMinRatio * kDelta) @@ -954,7 +974,7 @@ void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, Ali // // - debug part of FillResolutionHistoLocal - // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data - // called only for fDebugLevel > 4 + // called only for GetStreamLevel() > 4 // fill resolution trees // @@ -966,278 +986,276 @@ void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, Ali 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) ); - // - Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) ); - Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) ); - Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), + 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) ); + // + Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) ); + Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) ); + Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) ); - Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), + Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) ); - Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())), - TMath::Abs(anglez), TMath::Abs(cluster->GetMax()), - TMath::Sqrt(cluster0->GetSigmaY2()), 0 ); - Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())), - TMath::Abs(anglez), TMath::Abs(cluster->GetMax()), - TMath::Sqrt(cluster0->GetSigmaZ2()),0 ); - sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.); - sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.); - // - if (kDelta != 0){ - fitY0.AddPoint(&x, cluster->GetY(), sigmaY0); - fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0); - } - Double_t xxx[4]; - xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0; - xxx[1] = x; - xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0; - xxx[3] = x * x; - fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0); - fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ); - fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS); - fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0); - fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ); - fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS); - // - } // neigbouhood-loop - // - fitY0.Eval(); - fitZ0.Eval(); - fitY2.Eval(); - fitZ2.Eval(); - fitY2Q.Eval(); - fitZ2Q.Eval(); - fitY2S.Eval(); - fitZ2S.Eval(); - Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.); - Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.); - Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.); - Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.); - Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.); - Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.); - Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.); - Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.); - // - static TVectorD parY0(3); - static TMatrixD matY0(3, 3); - static TVectorD parZ0(3); - static TMatrixD matZ0(3, 3); - fitY0.GetParameters(parY0); - fitY0.GetCovarianceMatrix(matY0); - fitZ0.GetParameters(parZ0); - fitZ0.GetCovarianceMatrix(matZ0); - // - static TVectorD parY2(5); - static TMatrixD matY2(5,5); - static TVectorD parZ2(5); - static TMatrixD matZ2(5,5); - fitY2.GetParameters(parY2); - fitY2.GetCovarianceMatrix(matY2); - fitZ2.GetParameters(parZ2); - fitZ2.GetCovarianceMatrix(matZ2); - // - static TVectorD parY2Q(5); - static TMatrixD matY2Q(5,5); - static TVectorD parZ2Q(5); - static TMatrixD matZ2Q(5,5); - fitY2Q.GetParameters(parY2Q); - fitY2Q.GetCovarianceMatrix(matY2Q); - fitZ2Q.GetParameters(parZ2Q); - fitZ2Q.GetCovarianceMatrix(matZ2Q); - static TVectorD parY2S(5); - static TMatrixD matY2S(5,5); - static TVectorD parZ2S(5); - static TMatrixD matZ2S(5,5); - fitY2S.GetParameters(parY2S); - fitY2S.GetCovarianceMatrix(matY2S); - fitZ2S.GetParameters(parZ2S); - fitZ2S.GetCovarianceMatrix(matZ2S); - Float_t sigmaY0 = TMath::Sqrt(matY0(0,0)); - Float_t sigmaZ0 = TMath::Sqrt(matZ0(0,0)); - Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1)); - Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1)); - Float_t sigmaY2 = TMath::Sqrt(matY2(1,1)); - Float_t sigmaZ2 = TMath::Sqrt(matZ2(1,1)); - Float_t sigmaDY2 = TMath::Sqrt(matY2(3,3)); - Float_t sigmaDZ2 = TMath::Sqrt(matZ2(3,3)); - Float_t sigmaY2Q = TMath::Sqrt(matY2Q(1,1)); - Float_t sigmaZ2Q = TMath::Sqrt(matZ2Q(1,1)); - Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3)); - Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3)); - Float_t sigmaY2S = TMath::Sqrt(matY2S(1,1)); - Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1)); - Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3)); - Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3)); - - // Error parameters - Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley)); - Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez)); + Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster->GetMax()), + TMath::Sqrt(cluster0->GetSigmaY2()), 0 ); + Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster->GetMax()), + TMath::Sqrt(cluster0->GetSigmaZ2()),0 ); + sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.); + sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.); + // + if (kDelta != 0){ + fitY0.AddPoint(&x, cluster->GetY(), sigmaY0); + fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0); + } + Double_t xxx[4]; + xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0; + xxx[1] = x; + xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0; + xxx[3] = x * x; + fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0); + fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ); + fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS); + fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0); + fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ); + fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS); + // + } // neigbouhood-loop // - Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); - Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(anglez),TMath::Abs(cluster0->GetMax())); - Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); - Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(anglez),TMath::Abs(cluster0->GetMax())); - /////////////////////////////////////////////////////////////////////////////// -// // -// Class to analyse tracks for calibration // -// to be used as a component in AliTPCSelectorTracks // -// In the constructor you have to specify name and title // -// to get the Object out of a file. // -// The parameter 'clusterParam', a AliTPCClusterParam object // -// (needed for TPC cluster error and shape parameterization) // -// Normally you get this object out of the file 'TPCClusterParam.root' // -// In the parameter 'cuts' the cuts are specified, that decide // -// weather a track will be accepted for calibration or not. // -// // -// // -/////////////////////////////////////////////////////////////////////////////// - - // RMS parameters - Float_t meanRMSY = 0; - Float_t meanRMSZ = 0; - Int_t nclRMS = 0; - for (Int_t idelta = -2; idelta <= 2; idelta++){ - // loop over neighbourhood - if (idelta+irow < 0 || idelta+irow > 159) continue; -// if (idelta+irow>159) continue; - AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta); - if (!cluster) continue; - meanRMSY += TMath::Sqrt(cluster->GetSigmaY2()); - meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2()); - nclRMS++; - } - meanRMSY /= nclRMS; - meanRMSZ /= nclRMS; - - Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2()); - Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2()); - Float_t rmsYT = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); - Float_t rmsZT = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(anglez), TMath::Abs(cluster0->GetMax())); - Float_t rmsYT0 = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(angley)); - Float_t rmsZT0 = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + fitY0.Eval(); + fitZ0.Eval(); + fitY2.Eval(); + fitZ2.Eval(); + fitY2Q.Eval(); + fitZ2Q.Eval(); + fitY2S.Eval(); + fitZ2S.Eval(); + Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.); + Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.); + Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.); + Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.); + Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.); + Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.); + Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.); + Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.); + // + static TVectorD parY0(3); + static TMatrixD matY0(3, 3); + static TVectorD parZ0(3); + static TMatrixD matZ0(3, 3); + fitY0.GetParameters(parY0); + fitY0.GetCovarianceMatrix(matY0); + fitZ0.GetParameters(parZ0); + fitZ0.GetCovarianceMatrix(matZ0); + // + static TVectorD parY2(5); + static TMatrixD matY2(5,5); + static TVectorD parZ2(5); + static TMatrixD matZ2(5,5); + fitY2.GetParameters(parY2); + fitY2.GetCovarianceMatrix(matY2); + fitZ2.GetParameters(parZ2); + fitZ2.GetCovarianceMatrix(matZ2); + // + static TVectorD parY2Q(5); + static TMatrixD matY2Q(5,5); + static TVectorD parZ2Q(5); + static TMatrixD matZ2Q(5,5); + fitY2Q.GetParameters(parY2Q); + fitY2Q.GetCovarianceMatrix(matY2Q); + fitZ2Q.GetParameters(parZ2Q); + fitZ2Q.GetCovarianceMatrix(matZ2Q); + static TVectorD parY2S(5); + static TMatrixD matY2S(5,5); + static TVectorD parZ2S(5); + static TMatrixD matZ2S(5,5); + fitY2S.GetParameters(parY2S); + fitY2S.GetCovarianceMatrix(matY2S); + fitZ2S.GetParameters(parZ2S); + fitZ2S.GetCovarianceMatrix(matZ2S); + Float_t sigmaY0 = TMath::Sqrt(matY0(0,0)); + Float_t sigmaZ0 = TMath::Sqrt(matZ0(0,0)); + Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1)); + Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1)); + Float_t sigmaY2 = TMath::Sqrt(matY2(1,1)); + Float_t sigmaZ2 = TMath::Sqrt(matZ2(1,1)); + Float_t sigmaDY2 = TMath::Sqrt(matY2(3,3)); + Float_t sigmaDZ2 = TMath::Sqrt(matZ2(3,3)); + Float_t sigmaY2Q = TMath::Sqrt(matY2Q(1,1)); + Float_t sigmaZ2Q = TMath::Sqrt(matZ2Q(1,1)); + Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3)); + Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3)); + Float_t sigmaY2S = TMath::Sqrt(matY2S(1,1)); + Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1)); + Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3)); + Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3)); + + // Error parameters + Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley)); + Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez)); + // + Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); + Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(anglez),TMath::Abs(cluster0->GetMax())); + Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); + Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(anglez),TMath::Abs(cluster0->GetMax())); + + + // RMS parameters + Float_t meanRMSY = 0; + Float_t meanRMSZ = 0; + Int_t nclRMS = 0; + for (Int_t idelta = -2; idelta <= 2; idelta++){ + // loop over neighbourhood + if (idelta+irow < 0 || idelta+irow > 159) continue; + // if (idelta+irow>159) continue; + AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta); + if (!cluster) continue; + meanRMSY += TMath::Sqrt(cluster->GetSigmaY2()); + meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2()); + nclRMS++; + } + meanRMSY /= nclRMS; + meanRMSZ /= nclRMS; + + Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2()); + Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2()); + Float_t rmsYT = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(angley), TMath::Abs(cluster0->GetMax())); + Float_t rmsZT = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster0->GetMax())); + Float_t rmsYT0 = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(angley)); + Float_t rmsZT0 = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())), TMath::Abs(anglez)); - Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(anglez), TMath::Abs(cluster0->GetMax())); - Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(anglez), TMath::Abs(cluster0->GetMax())); - Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()), - rmsY,meanRMSY); - Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), - TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()), - rmsZ,meanRMSZ); - - // cluster debug - (*fDebugStream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly - "Sector="<GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster0->GetMax())); + Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster0->GetMax())); + Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()), + rmsY,meanRMSY); + Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())), + TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()), + rmsZ,meanRMSZ); + // + // cluster debug + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + (*cstream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly + "run="< 6) Info("Draw", "Drawing an exemplaric picture."); + if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture."); SetStyle(); Double_t min = 0; Double_t max = 0; @@ -1345,7 +1363,7 @@ void AliTPCcalibTracks::Draw(Option_t* 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' @@ -1353,7 +1371,7 @@ 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); + if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName); MakeAmpPlots(stat, pathName); MakeDeltaPlots(pathName); FitResolutionNew(pathName); @@ -1364,7 +1382,7 @@ void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ } -void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ +void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){ // // creates several plots: // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps @@ -1391,7 +1409,7 @@ void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ allAmpHisOROC->SetTitle("Amp all OROCs"); ps = new TPostScript("fArrayAmp.ps", 112); - if (fDebugLevel > -1) cout << "creating fArrayAmp.ps..." << endl; + if (GetDebugLevel() > -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(); @@ -1415,7 +1433,7 @@ void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ Double_t min = 0; Double_t max = 0; ps = new TPostScript("fArrayAmpRow.ps", 112); - if (fDebugLevel > -1) cout << "creating fArrayAmpRow.ps..." << endl; + if (GetDebugLevel() > -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; @@ -1434,7 +1452,7 @@ void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ } -void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){ +void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){ // // creates several plots: // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit @@ -1451,7 +1469,7 @@ void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){ Double_t max = 0; ps = new TPostScript("DeltaYZ.ps", 112); - if (fDebugLevel > -1) cout << "creating DeltaYZ.ps..." << endl; + if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl; min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20; max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20; fDeltaY->SetAxisRange(min, max); @@ -1471,7 +1489,7 @@ void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){ } -void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){ +void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const 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 @@ -1484,7 +1502,7 @@ void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* 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; + if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl; TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone()); TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone()); chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC"); @@ -1513,7 +1531,7 @@ void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){ } -void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){ +void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){ // // creates charge vs. driftlength plots, one TProfile for each ROC // under development.... @@ -1528,7 +1546,7 @@ void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){ c1->Divide(0,3); TPostScript *ps; ps = new TPostScript("chargeVsDriftlength.ps", 111); - if (fDebugLevel > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl; + if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl; TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone()); TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone()); @@ -1575,7 +1593,7 @@ void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){ -void AliTPCcalibTracks::FitResolutionNew(char* pathName){ +void AliTPCcalibTracks::FitResolutionNew(const char* pathName){ // // calculates different resulution fits in Y and Z direction // the histograms are written to 'ResolutionYZ.ps' @@ -1589,7 +1607,7 @@ void AliTPCcalibTracks::FitResolutionNew(char* pathName){ TCanvas c; c.Divide(2,1); - if (fDebugLevel > -1) cout << "creating ResolutionYZ.ps..." << endl; + if (GetDebugLevel() > -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); @@ -1698,7 +1716,7 @@ void AliTPCcalibTracks::FitResolutionNew(char* pathName){ } -void AliTPCcalibTracks::FitRMSNew(char* pathName){ +void AliTPCcalibTracks::FitRMSNew(const char* pathName){ // // calculates different resulution-rms fits in Y and Z direction // the histograms are written to 'RMS_YZ.ps' @@ -1712,7 +1730,7 @@ void AliTPCcalibTracks::FitRMSNew(char* 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; + if (GetDebugLevel() > -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); @@ -1821,7 +1839,7 @@ void AliTPCcalibTracks::FitRMSNew(char* pathName){ } -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' @@ -1849,8 +1867,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); @@ -1884,7 +1902,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; //-------------------------------------------------------------------------------------------- @@ -1904,7 +1922,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; @@ -2098,7 +2116,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); @@ -2121,8 +2139,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(".."); } @@ -2130,384 +2148,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 @@ -2515,12 +2155,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 @@ -2550,9 +2190,9 @@ 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; @@ -2575,15 +2215,17 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { 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; } // merge data members - if (fDebugLevel > 0) cout << "histogram's merge-functins are called... " << endl; + if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl; fDeltaY->Merge(deltaYList); fDeltaZ->Merge(deltaZList); fHclus->Merge(hclusList); @@ -2597,7 +2239,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { TList* histList = 0; TIterator *objListIterator = 0; - if (fDebugLevel > 0) cout << "merging fArrayAmpRows..." << endl; + if (GetDebugLevel() > 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(); @@ -2613,12 +2255,12 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete objListIterator; } - if (fDebugLevel > 0) cout << "merging fArrayAmps..." << endl; + if (GetDebugLevel() > 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(); + TIterator *cobjListIterator = arrayAmpList->MakeIterator(); histList = new TList; - while (( objarray = (TObjArray*)objListIterator->Next() )) { + while (( objarray = (TObjArray*)cobjListIterator->Next() )) { // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F if (!objarray) continue; hist = (TH1F*)(objarray->At(i)); @@ -2626,10 +2268,10 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { } ((TH1F*)(fArrayAmp->At(i)))->Merge(histList); delete histList; - delete objListIterator; + delete cobjListIterator; } - 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(); @@ -2645,7 +2287,7 @@ 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(); @@ -2661,7 +2303,7 @@ 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(); @@ -2677,7 +2319,7 @@ 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(); @@ -2693,7 +2335,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete objListIterator; } - if (fDebugLevel > 0) cout << "merging fArrayChargeVsDriftlength..." << endl; + if (GetDebugLevel() > 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(); @@ -2709,7 +2351,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { delete objListIterator; } - if (fDebugLevel > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl; + if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl; // merge fcalPadRegionChargeVsDriftlength AliTPCCalPadRegion *cpr = 0x0; @@ -2744,7 +2386,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { - 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(); @@ -2817,12 +2459,9 @@ 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; } @@ -2882,5 +2521,124 @@ AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClu } +void AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){ + // + // Make position corrections + // for the moment Only using debug streamer + // chainres - debug tree + // param - parameters to be updated + // maxPoints - maximal number of points using for fit + // verbose - print info flag + // + // Current implementation - only using debug streamers + // + + /* + //Defaults + Int_t maxPoints=100000; + */ + /* + Usage: + //0. Load libraries + gSystem->Load("libANALYSIS"); + gSystem->Load("libSTAT"); + gSystem->Load("libTPCcalib"); + + + //1. Load Parameters to be modified: + //e.g: + + AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); + AliCDBManager::Instance()->SetRun(0) + AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam(); + + //2. Load chain from debug streamers + // + //e.g + gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); + gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+") + AliXRDPROOFtoolkit tool; + TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200); + chainres->Lookup(); + //3. Do fits and store results + // + AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ; + TFile f("paramout.root","recreate"); + param->Write("clusterParam"); + f.Close(); + //4. Verification + TFile f2("paramout.root"); + AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam"); + param2->SetInstance(param2); + chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line + + */ + + + TStatToolkit toolkit; + Double_t chi2; + TVectorD fitParamY0; + TVectorD fitParamY1; + TVectorD fitParamZ0; + TVectorD fitParamZ1; + TMatrixD covMatrix; + Int_t npoints; + + chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)"); + chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)"); + chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)"); + chainres->SetAlias("st","(sin(dt)-dt)"); + // + chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))"); + // + // + // + TCut cutA("1"); + TString fstringY=""; + // + fstringY+="(dp)++"; //1 + fstringY+="(dp)*di++"; //2 + fstringY+="(sp)++"; //3 + fstringY+="(sp)*di++"; //4 + TString fstringZ=""; + fstringZ+="(dt)++"; //1 + fstringZ+="(dt)*di++"; //2 + fstringZ+="(st)++"; //3 + fstringZ+="(st)*di++"; //4 + // + // Z corrections + // + TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints); + printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); + param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone(); + // + TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints); + // + printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); + param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone(); + param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone(); + // + // Y corrections + // + TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints); + printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); + param->PosYcor(0) = (TVectorD*) fitParamY0.Clone(); + + + TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints); + // + printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); + param->PosYcor(1) = (TVectorD*) fitParamY1.Clone(); + param->PosYcor(2) = (TVectorD*) fitParamY1.Clone(); + // + // + // + chainres->SetAlias("fitZ0",strZ0->Data()); + chainres->SetAlias("fitZ1",strZ1->Data()); + chainres->SetAlias("fitY0",strY0->Data()); + chainres->SetAlias("fitY1",strY1->Data()); + // chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000); +} +