From 10757ee9226dad26cb3b3f26602615a3cbc31767 Mon Sep 17 00:00:00 2001 From: marian Date: Mon, 3 Dec 2007 16:27:21 +0000 Subject: [PATCH] Moving the TPC calibration using tracks from separate directory to the TPC directory The TPC calibration will be compiled by default aliroot compilation (Marian) --- TPC/AliTPCCalPadRegion.cxx | 117 ++ TPC/AliTPCCalPadRegion.h | 48 + TPC/AliTPCFitPad.cxx | 143 ++ TPC/AliTPCFitPad.h | 40 + TPC/AliTPCSelectorESD.cxx | 431 +++++ TPC/AliTPCSelectorESD.h | 67 + TPC/AliTPCSelectorTracks.cxx | 228 +++ TPC/AliTPCSelectorTracks.h | 43 + TPC/AliTPCcalibTracks.cxx | 2841 +++++++++++++++++++++++++++++++++ TPC/AliTPCcalibTracks.h | 167 ++ TPC/AliTPCcalibTracksCuts.cxx | 148 ++ TPC/AliTPCcalibTracksGain.cxx | 983 ++++++++++++ TPC/AliTPCcalibTracksGain.h | 110 ++ TPC/AliTPCcalibV0.cxx | 811 ++++++++++ TPC/AliTPCcalibV0.h | 66 + TPC/TPCcalibLinkDef.h | 28 + TPC/libTPCcalib.pkg | 14 + 17 files changed, 6285 insertions(+) create mode 100644 TPC/AliTPCCalPadRegion.cxx create mode 100644 TPC/AliTPCCalPadRegion.h create mode 100644 TPC/AliTPCFitPad.cxx create mode 100644 TPC/AliTPCFitPad.h create mode 100644 TPC/AliTPCSelectorESD.cxx create mode 100644 TPC/AliTPCSelectorESD.h create mode 100644 TPC/AliTPCSelectorTracks.cxx create mode 100644 TPC/AliTPCSelectorTracks.h create mode 100644 TPC/AliTPCcalibTracks.cxx create mode 100644 TPC/AliTPCcalibTracks.h create mode 100644 TPC/AliTPCcalibTracksCuts.cxx create mode 100644 TPC/AliTPCcalibTracksGain.cxx create mode 100644 TPC/AliTPCcalibTracksGain.h create mode 100644 TPC/AliTPCcalibV0.cxx create mode 100644 TPC/AliTPCcalibV0.h create mode 100644 TPC/TPCcalibLinkDef.h create mode 100644 TPC/libTPCcalib.pkg diff --git a/TPC/AliTPCCalPadRegion.cxx b/TPC/AliTPCCalPadRegion.cxx new file mode 100644 index 00000000000..7a8e67284a4 --- /dev/null +++ b/TPC/AliTPCCalPadRegion.cxx @@ -0,0 +1,117 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//////////////////////////////////////////////////////////////////////////// +// +// === Class for properties specific to pad regions === +// +// Each segment of the TPC (i.e. IROC and corresponding OROC) consists +// of three different pad sizes (short, medium, long). This class +// is useful for scenarios, where it is appropriate to describe +// some behaviour per pad size region. It provides an easy interface +// for getting and setting arbitrary objects for each pad size region. +// There is no need that this object is of the same type for each +// pad size region (though it probably will be in most of the cases), +// nor that it is set at all (e.g. when no data for this region +// exists and such an object is not needed). +// +// An example that makes usage of this class is the AliTPCFitPad class +// which stores TLinearFitter objects for each pad region. +// +//////////////////////////////////////////////////////////////////////////// + +#include "AliTPCCalPadRegion.h" +#include "AliTPCROC.h" + +ClassImp(AliTPCCalPadRegion) + +AliTPCCalPadRegion::AliTPCCalPadRegion(): + TNamed(), + fObjects(0) +{ + // + // Default constructor. + // +} + +AliTPCCalPadRegion::AliTPCCalPadRegion(const char *name, const char *title) : + TNamed(name, title), + fObjects(0) +{ + // + // Constructor. + // + + fObjects = new TObjArray(fgkNSegments * fgkNPadTypes); +} + +AliTPCCalPadRegion::AliTPCCalPadRegion(const AliTPCCalPadRegion& obj) : + TNamed(obj) +{ + // + // Copy constructor. + // + + fObjects = new TObjArray(*(obj.fObjects)); +} + +AliTPCCalPadRegion& AliTPCCalPadRegion::operator=(const AliTPCCalPadRegion& rhs) { + // + // Assignment operator. + // + + if (this != &rhs) { + TNamed::operator=(rhs); + fObjects = new TObjArray(*(rhs.fObjects)); + } + return *this; +} + +void AliTPCCalPadRegion::GetPadRegionCenterLocal(UInt_t padType, Double_t* xy) { + // + // Return the center of the pad size region in local + // coordinates as an Double_t array xy of length 2. + // + + Float_t centerPad[3] = {0}; + AliTPCROC* tpcROC = AliTPCROC::Instance(); + + Int_t IOROC = (padType == 0) ? 0 : tpcROC->GetNInnerSector(); + //tpcROC->GetPositionLocal(IOROC, tpcROC->GetNRows(IOROC)/2, tpcROC->GetNPads(IOROC, tpcROC->GetNRows(IOROC)/2)/2, centerPad); // use this instead of the switch statement if you want to calculate the center of the ROC and not the center of the regions with the same pad size + switch (padType) { + case 0: // short pads + tpcROC->GetPositionLocal(IOROC, tpcROC->GetNRows(IOROC)/2, tpcROC->GetNPads(IOROC, tpcROC->GetNRows(IOROC)/2)/2, centerPad); + break; + case 1: // medium pads + tpcROC->GetPositionLocal(IOROC, 64/2, tpcROC->GetNPads(IOROC, 64/2)/2, centerPad); + break; + case 2: // long pads + tpcROC->GetPositionLocal(IOROC, 64+32/2, tpcROC->GetNPads(IOROC, 64+32/2)/2, centerPad); + break; + } + + xy[0] = centerPad[0]; + xy[1] = centerPad[1]; +} + +/*UInt_t AliTPCCalPadRegion::GetStartRow(UInt_t padType) { + // + // Returns the index of the + // +} + +UInt_t AliTPCCalPadRegion::GetEndRow(UInt_t padType) { + +}*/ diff --git a/TPC/AliTPCCalPadRegion.h b/TPC/AliTPCCalPadRegion.h new file mode 100644 index 00000000000..25a6f03e783 --- /dev/null +++ b/TPC/AliTPCCalPadRegion.h @@ -0,0 +1,48 @@ +#ifndef ALITPCCALPADREGION_H +#define ALITPCCALPADREGION_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +#include +#include +#include +#include + +class AliTPCCalPadRegion: public TNamed { +public: + AliTPCCalPadRegion(); + AliTPCCalPadRegion(const AliTPCCalPadRegion& obj); + AliTPCCalPadRegion(const char *name, const char *title); + //AliTPCCalPadRegion(const TString &name, const TString &title) : TNamed(name, title) { } + virtual ~AliTPCCalPadRegion() { delete fObjects; } + AliTPCCalPadRegion& operator=(const AliTPCCalPadRegion& rhs); + + virtual TObject* GetObject(UInt_t segment, UInt_t padType) + { return BoundsOk("GetObject", segment, padType) ? fObjects->At(segment+fgkNSegments*padType) : 0x0; } + virtual void SetObject(TObject* obj, UInt_t segment, UInt_t padType) + { if (BoundsOk("SetObject", segment, padType)) fObjects->AddAt(obj, segment+fgkNSegments*padType); } + virtual void Delete(Option_t* option = "") { if (fObjects) fObjects->Delete(); } + virtual TIterator* MakeIterator(Bool_t direction = kIterForward) const { return fObjects->MakeIterator(direction); } + static UInt_t GetNSegments() { return fgkNSegments; } + static UInt_t GetNPadTypes() { return fgkNPadTypes; } + static void GetPadRegionCenterLocal(UInt_t padType, Double_t* xy); +// static UInt_t GetStartRow(UInt_t padType); +// static UInt_t GetEndRow(UInt_t padType); + +protected: + virtual Bool_t BoundsOk(const char* where, UInt_t segment, UInt_t padType) const + { return (segment >= fgkNSegments || padType >= fgkNPadTypes) ? OutOfBoundsError(where, segment, padType) : kTRUE; } + virtual Bool_t OutOfBoundsError(const char* where, UInt_t segment, UInt_t padType) const + { Error(where, "Index out of bounds (trying to access segment %d, pad type %d).", segment, padType); return kFALSE; } + + TObjArray* fObjects; // array containing an object for each pad region + + static const UInt_t fgkNSegments = 36; // number of TPC sectors, 0-17: A side, 18-35: C side (IROC and OROC are treated as one sector) + static const UInt_t fgkNPadTypes = 3; // number of pad types, 0: short pads, 1: medium pads, 2: long pads + + ClassDef(AliTPCCalPadRegion, 1) +}; + + +#endif diff --git a/TPC/AliTPCFitPad.cxx b/TPC/AliTPCFitPad.cxx new file mode 100644 index 00000000000..a9d1022e1fa --- /dev/null +++ b/TPC/AliTPCFitPad.cxx @@ -0,0 +1,143 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//////////////////////////////////////////////////////////////////////////// +// +// === Class for fitting properties specific to pad regions === +// +// For each pad size region a separate TLinearFitter object is assigned. +// Commonly used functions such as getting the center of a pad size +// region or visualization functions are provided. Also, choosing the +// segment and pad type is made easy due to different methods that +// can calculate these informations from other coordinates. +// +//////////////////////////////////////////////////////////////////////////// + +#include "AliTPCFitPad.h" + +ClassImp(AliTPCFitPad) + +AliTPCFitPad::~AliTPCFitPad() { + // + // Destructor. + // + + Delete(); +} + +AliTPCFitPad::AliTPCFitPad(Int_t ndim, const char* formula, Option_t* opt) : + AliTPCCalPadRegion("", ""), + fNdim(ndim), + fFormula(formula), + fOpt(opt) +{ + // + // Constructor. The parameters are used for generating new TLinearFitter + // objects and are described in its documentation. + // +} + +AliTPCFitPad& AliTPCFitPad::operator=(const AliTPCFitPad& rhs) { + // + // Assignment operator. + // + + if (this != &rhs) { + AliTPCCalPadRegion::operator=(rhs); + fNdim = rhs.fNdim; + fFormula = rhs.fFormula; + fOpt = rhs.fOpt; + } + return *this; +} + +void AliTPCFitPad::Add(AliTPCFitPad* fit) { + // + // Adds another AliTPCFitPad object to this object. The formula should be the + // same, though it won't be checked! + // + + for (UInt_t iSegment = 0; iSegment < GetNSegments(); iSegment++) { + for (UInt_t iPadType = 0; iPadType < GetNPadTypes(); iPadType++) { + TLinearFitter* fitter = fit->GetFitterSimple(iSegment, iPadType); + // parameter workaround == kTRUE because it is not possible to add another + // TLinearFitter object to a "virgin" one. Thus a dummy data point is added + // and cleared again immediately afterwards. Due to a broken TLinearFitter + // copy constructor this is a necessary workaround. + if (fitter) { + cout << "TLinearFitter::Add called for " << iSegment << ", " << iPadType << endl; + GetFitter(iSegment, iPadType, kTRUE)->Add(fitter); + } + } + } +} + +TLinearFitter* AliTPCFitPad::GetFitterSimple(UInt_t segment, UInt_t padType) { + // + // This method returns the fitter corresponding to segment and pad type. + // In contrast to GetFitter() no fitter will be created, if it does + // not exist, but a null pointer is returned. + // + + return (TLinearFitter*)(GetObject(segment, padType)); +} + +TLinearFitter* AliTPCFitPad::GetFitter(UInt_t segment, UInt_t padType, Bool_t workaround) { + // + // This method returns the fitter corresponding + // to segment and pad type. + // If the fitter doesn't exist yet, it will be created on the fly + // according to the parameters passed to the constructor. + // + // The workaround parameter should always be kFALSE. (It is only used by the Add method.) + // + + TLinearFitter* fitter = GetFitterSimple(segment, padType); + if (fitter == 0) { + SetObject(new TLinearFitter(fNdim, fFormula, fOpt), segment, padType); + fitter = (TLinearFitter*)(GetObject(segment, padType)); + if (workaround) { + Double_t x[1000]; + for (Int_t i = 0; i < fNdim; i++) x[i] = 3.141592; + fitter->AddPoint(x, 31.41592); + fitter->ClearPoints(); + //cout << "workaround called for " << segment << ", " << padType << endl; + } + } + return fitter; +} + +Int_t AliTPCFitPad::Evaluate(Bool_t robust, Double_t frac) { + // + // Evaluates all fitters. Returns 0 if successful, 1 in case of an error. + // If the robust option is set to kTRUE a robust fit is performed with frac as + // the minimal fraction of good points (see TLinearFitter::EvalRobust for details). + // Beware: Robust fitting is much slower! + // + + Int_t returnCode = 0; + for (UInt_t iSegment = 0; iSegment < GetNSegments(); iSegment++) { + for (UInt_t iPadType = 0; iPadType < GetNPadTypes(); iPadType++) { + if (TLinearFitter* fitter = GetFitterSimple(iSegment, iPadType)) { + Int_t status = robust ? fitter->EvalRobust(frac) : fitter->Eval(); + if (status != 0) { + returnCode = 1; + Error("Evaluate", "Error in evaluation of fitter in segment %d, pad region %d", iSegment, iPadType); + } + } + } + } + return returnCode; +} diff --git a/TPC/AliTPCFitPad.h b/TPC/AliTPCFitPad.h new file mode 100644 index 00000000000..06d51e3244d --- /dev/null +++ b/TPC/AliTPCFitPad.h @@ -0,0 +1,40 @@ +#ifndef ALITPCFITPAD_H +#define ALITPCFITPAD_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +#include +#include "AliTPCCalPadRegion.h" +#include +#include + +using namespace std; + +class TString; + +class AliTPCFitPad: public AliTPCCalPadRegion { +public: + AliTPCFitPad() : AliTPCCalPadRegion() { } + AliTPCFitPad(const AliTPCFitPad& obj) : AliTPCCalPadRegion(obj), fNdim(obj.fNdim), fFormula(obj.fFormula), fOpt(obj.fOpt) { } + AliTPCFitPad(Int_t ndim, const char* formula, Option_t* opt = "D"); + AliTPCFitPad& operator=(const AliTPCFitPad& rhs); + //AliTPCFitPad(const char *name, const char *title) : AliTPCCalPadRegion(name, title) { } + //AliTPCFitPad(const TString &name, const TString &title) : AliTPCCalPadRegion(name, title) { } + virtual ~AliTPCFitPad(); + + void Add(AliTPCFitPad* fit); + TLinearFitter* GetFitter(UInt_t segment, UInt_t padType, Bool_t workaround = kFALSE); + TLinearFitter* GetFitterSimple(UInt_t segment, UInt_t padType); + Int_t Evaluate(Bool_t robust = kFALSE, Double_t frac = -1.); + +protected: + Int_t fNdim; // used for generating new TLinearFitter objects + TString fFormula; // used for generating new TLinearFitter objects + TString fOpt; // used for generating new TLinearFitter objects + + ClassDef(AliTPCFitPad, 1) +}; + + +#endif diff --git a/TPC/AliTPCSelectorESD.cxx b/TPC/AliTPCSelectorESD.cxx new file mode 100644 index 00000000000..5bebfb6a2f2 --- /dev/null +++ b/TPC/AliTPCSelectorESD.cxx @@ -0,0 +1,431 @@ +// The class definition in esdClus.h has been generated automatically +// by the ROOT utility TTree::MakeSelector(). This class is derived +// from the ROOT class TSelector. For more information on the TSelector +// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual. + +// The following methods are defined in this file: +// Begin(): called everytime a loop on the tree starts, +// a convenient place to create your histograms. +// SlaveBegin(): called after Begin(), when on PROOF called only on the +// slave servers. +// Process(): called for each event, in this function you decide what +// to read and fill your histograms. +// SlaveTerminate: called at the end of the loop on the tree, when on PROOF +// called only on the slave servers. +// Terminate(): called at the end of the loop on the tree, +// a convenient place to draw/fit your histograms. +// +// To use this file, try the following session on your Tree T: +// +// Root > T->Process("esdClus.C") +// Root > T->Process("esdClus.C","some options") +// Root > T->Process("esdClus.C+") +// +// Modification log: +// 05/11/2006 HH Correct for large pads (outer sectors) in amplitude plots + +#include "TSystem.h" +#include +#include +#include "TCint.h" +#include "TH1I.h" +#include "TTimeStamp.h" +#include "TProof.h" +#include "TTree.h" +// +#include "AliTracker.h" +#include "AliMagF.h" +// +#include "AliESDEvent.h" // new container +#include "AliESD.h" +#include "AliESDtrack.h" +#include "AliESDfriend.h" +#include "AliESDfriendTrack.h" +#include "AliTPCseed.h" +#include "AliTPCclusterMI.h" +// +#include "AliSysInfo.h" +#include "AliTPCSelectorESD.h" + + + + + +AliTPCSelectorESD::AliTPCSelectorESD(TTree *tree) : + TSelector(), + fChain(0), + fESDevent(0), + fESD(0), + fESDfriend(0), + fFileNo(0), + fSysWatch(0), // system info + fFileWatch(0), // file watch - write the status of the analyzed files + fDebugLevel(0) + { + G__SetCatchException(0); + // + // + //if (somthing){ + fSysWatch = new fstream("syswatch2.log", ios_base::out|ios_base::trunc); + fFileWatch = new fstream("filewatch.log", ios_base::out|ios_base::trunc); + if (gProof) fDebugLevel = gProof->GetLogLevel(); + if (tree) fChain=tree; + } + + +void AliTPCSelectorESD::Begin(TTree * /*tree*/) +{ + // The Begin() function is called at the start of the query. + // When running with PROOF Begin() is only called on the client. + // The tree argument is deprecated (on PROOF 0 is passed). + + TString option = GetOption(); + +} + +void AliTPCSelectorESD::SlaveBegin(TTree * tree) +{ + // The SlaveBegin() function is called after the Begin() function. + // When running with PROOF SlaveBegin() is called on each slave server. + // The tree argument is deprecated (on PROOF 0 is passed). + if (tree) fChain = tree; + Init(tree); + // + fNtracks = new TH1I("ntracks","Number of tracks",100,0,400); + fNtracksFriend = new TH1I("ntracksF","Number of friend tracks",100,0,400); + fNClusters = new TH1I("ncluster","Number of clusters",100,0,200); + fOutput->AddLast(fNtracks); + fOutput->AddLast(fNtracksFriend); + fOutput->AddLast(fNClusters); + + +} + +void AliTPCSelectorESD::CleanESD(){ + // + Bool_t isNew = fESDevent!=0; + if (isNew) return; + + + if (fESD!=0){ + delete fESD; + fESD = 0; + } + if (fESDevent!=0){ + delete fESDevent; + fESDevent=0; + } + if (fESDfriend){ + delete fESDfriend; + fESDfriend =0; + } +} + +Bool_t AliTPCSelectorESD::Process(Long64_t entry) +{ + // The Process() function is called for each entry in the tree (or possibly + // keyed object in the case of PROOF) to be processed. The entry argument + // specifies which entry in the currently loaded tree is to be processed. + // It can be passed to either AliTPCSelectorESD::GetEntry() or TBranch::GetEntry() + // to read either all or the required parts of the data. When processing + // keyed objects with PROOF, the object is already loaded and is available + // via the fObject pointer. + // + // This function should contain the "body" of the analysis. It can contain + // simple or elaborate selection criteria, run algorithms on the data + // of the event and typically fill histograms. + // + // The processing can be stopped by calling Abort(). + // + // Use fStatus to set the return value of TTree::Process(). + // + // The return value is currently not used. + + Int_t status = ProcessIn(entry); + if (fFileWatch) { + (*fFileWatch) << "__" << status ; + } + return (status==0); +} + + + + +Int_t AliTPCSelectorESD::ReadEvent(Long64_t entry){ + // + // + // + + + if (!fChain) return -1; + if (!fChain->GetTree()) return -1; + try { + fChain->GetTree()->GetEntry(entry); + } catch (std::bad_alloc) { + printf("Pica vyjebana pojebany skurveny kokot piciak\n"); + // fESD =0; + //fESDfriend = 0; + //fESDevent=0; + return -1; + } + // + // Info("Procces","0"); + if (!fESD && !fESDevent) { + //fESD = 0; + //fESDfriend = 0; + //CleanESD(); + return -2; + } + Int_t ntracks = (fESD) ? fESD->GetNumberOfTracks() : fESDevent->GetNumberOfTracks(); + + if (fESDevent){ + fESDevent->SetESDfriend(fESDfriend); + } + + + fNtracks->Fill(ntracks); + Info("Procces", Form("entry\t%d: Ntracks = %d",entry, ntracks)); + + if (!fESDfriend || fESDfriend->GetNumberOfTracks() != ntracks) { + try { + delete fESD; + } + catch (std::bad_alloc) { + printf("Pica vyjebana pojebany skurveny kokot piciak\n"); + fESD =0; + return -1; + } + //fESD = 0; + //fESDfriend = 0; + //fESD = 0; + // CleanESD(); + if (fESDfriend) fNtracksFriend->Fill(fESDfriend->GetNumberOfTracks()); + Info("Procces","2: PROBLEM"); + return -3; + } + if (fESD) fESD->SetESDfriend(fESDfriend); + return 0; +} + + +Int_t AliTPCSelectorESD::ProcessIn(Long64_t entry) +{ + // + // User part of proccess method + // + + // + // USER code to go here + // + Int_t status = ReadEvent(entry); + if (status<0) return status; + Int_t ntracks = (fESD) ? fESD->GetNumberOfTracks() : fESDevent->GetNumberOfTracks(); + // + AliTPCseed *seed; + AliTPCclusterMI cl; + Info("Procces", Form("entry\t%d: NtracksF = %d",entry,fESDfriend->GetNumberOfTracks() )); + + for (Int_t tr = 0; tr < ntracks; tr++){ + AliESDtrack *esdTrack = fESD ? (AliESDtrack*) fESD->GetTrack(tr): (AliESDtrack*) fESDevent->GetTrack(tr); + AliESDfriendTrack *friendtrack = (AliESDfriendTrack*) esdTrack->GetFriendTrack(); + seed = 0; + TObject *cobject=0; + for (Int_t i=0;;i++){ + cobject = friendtrack->GetCalibObject(i); + if (!cobject) break; + seed = dynamic_cast(cobject); + if (seed) break; + } + // + //Info("Process",Form("Proccessing track%d\n",tr)); + if (seed) { + fNClusters->Fill(seed->GetNumberOfClusters()); + // + // + } + } + CleanESD(); + return 0; + +} + + +void AliTPCSelectorESD::SlaveTerminate() +{ + // The SlaveTerminate() function is called after all entries or objects + // have been processed. When running with PROOF SlaveTerminate() is called + // on each slave server. + printf ("SlaveTerminate.. \n"); + +} + +void AliTPCSelectorESD::Terminate() +{ + // The Terminate() function is the last function to be called during + // a query. It always runs on the client, it can be used to present + // the results graphically or save the results to file. + + printf ("Terminate... \n"); + if (!fOutput) return; + TFile file("Output.root","recreate"); + printf("fOutput contains the following: \n"); + fOutput->Print(); + printf("Trying to write the file 'Output.root'... \n"); + fOutput->Write(); + file.Close(); + + +} + + + +void AliTPCSelectorESD::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normaly not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set branch addresses and branch pointers + static Int_t counter=0; + printf(Form("\nAliTPCSelectorESD::Init Accesing%d time\n",counter)); + counter++; + if (!tree) return; + fChain = tree; + //if (counter>1) return; + tree->SetBranchStatus("*",1); + // + // New AliESDevent format + // + if (!fChain->GetBranch("ESD")){ + // + // + // + if (fESDevent) delete fESDevent; + fESDevent = new AliESDEvent(); + fESDevent->ReadFromTree(tree); // Attach the branch with ESD friends + fESDfriend = (AliESDfriend*)fESDevent->FindListObject("AliESDfriend"); + tree->SetBranchAddress("ESDfriend.",&fESDfriend); + return; + } + // + // if old format + // + + + + // fChain->SetMakeClass(1); + fChain->SetBranchAddress("ESD",&fESD); + Info("Init","Enter"); + Bool_t isOK=kFALSE; + if (fChain->GetBranch("ESDfriend")) { + fChain->SetBranchAddress("ESDfriend",&fESDfriend); + Info("Init","V0-ESDfriend."); + isOK=kTRUE; + } + if (fChain->GetBranch("ESDfriend.")){ + Info("Init","V1-ESDfriend."); + fChain->SetBranchAddress("ESDfriend.",&fESDfriend); + isOK=kTRUE; + } + if (isOK) return; + + // + // Try to solve problem + // + + Info("Init","Problem"); + if (tree->GetBranch("ESD")){ + Info("InitTree",tree->GetBranch("ESD")->GetFile()->GetName()); + char fname[1000]; + sprintf(fname,"%s/AliESDfriends.root",gSystem->DirName(tree->GetBranch("ESD")->GetFile()->GetName())); + Info("InitFile",fname); + if (tree->AddFriend("esdFriendTree",fname)){ + Info("InitFileOK",fname); + if (fChain->GetBranch("ESDfriend")) { + fChain->SetBranchAddress("ESDfriend",&fESDfriend); + Info("Init","V0-ESDfriend."); + isOK=kTRUE; + } + if (fChain->GetBranch("ESDfriend.")){ + Info("Init","V1-ESDfriend."); + fChain->SetBranchAddress("ESDfriend.",&fESDfriend); + isOK=kTRUE; + } + } + } +} + + + +Bool_t AliTPCSelectorESD::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normaly not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + ++fFileNo; + const char * fname = "UNKNOWN"; + const char * hname = gSystem->HostName(); + if (!fChain) return kFALSE; + if (fChain->GetCurrentFile()){ + fname = fChain->GetCurrentFile()->GetName(); + } + Info("Notify",Form("Host %s processing file no %d %s\n",hname,fFileNo,fname)); + + // + // Print statistic to log file + // +// if (fname) { +// (*fFileWatch) << endl; +// (*fFileWatch) << hname <<"\t" +// << fname <<"\t"; +// } + DumpSysInfo(-1); + + return kTRUE; +} + + +void AliTPCSelectorESD::DumpSysInfo(Int_t entry){ + // + // dump system info to log file + // entry - entry number in the chain + // + const char * fname = "UNKNOWN"; + const char * hname = gSystem->HostName(); + if (fChain->GetCurrentFile()){ + fname = fChain->GetCurrentFile()->GetName(); + } + // // + if (fSysWatch){ + TTimeStamp stamp; + CpuInfo_t cpuInfo; + MemInfo_t memInfo; + ProcInfo_t procInfo; + + gSystem->GetCpuInfo(&cpuInfo, 1000); + gSystem->GetMemInfo(&memInfo); + gSystem->GetProcInfo(&procInfo); + + (*fSysWatch) << hname <<"\t" // hname - hostname + << fname <<"\t" // fname - filename + << entry <<"\t" // entry - entry number + << stamp.GetSec()<<"\t" // time - time stamp in seconds + << memInfo.fMemUsed<<"\t" // + << memInfo.fSwapUsed<<"\t" // + << procInfo.fMemResident<<"\t" // + << procInfo.fMemVirtual<<"\t" // + << cpuInfo.fUser <<"\t" // + << cpuInfo.fSys <<"\t" // + << procInfo.fCpuUser<<"\t" // + << procInfo.fCpuSys<<"\t" // + << endl; + } + AliSysInfo::AddStamp(fname); +} diff --git a/TPC/AliTPCSelectorESD.h b/TPC/AliTPCSelectorESD.h new file mode 100644 index 00000000000..c8955684d70 --- /dev/null +++ b/TPC/AliTPCSelectorESD.h @@ -0,0 +1,67 @@ + +#ifndef AliTPCSelectorESD_h +#define AliTPCSelectorESD_h + +#include +#include +using namespace std; +#include + +#include +#include +#include + +class AliESDEvent; +class AliESD; +class AliESDfriend; +class TH1I; + + +class AliTPCSelectorESD : public TSelector { +public : + AliTPCSelectorESD(TTree *tree=0); + virtual ~AliTPCSelectorESD() { /*delete fESD; delete fESDfriend;*/ } + virtual Int_t Version() const { return 1; } + virtual void Begin(TTree *tree); + virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); + virtual Bool_t Notify(); + virtual Bool_t Process(Long64_t entry); + virtual Int_t ReadEvent(Long64_t entry); + virtual Int_t ProcessIn(Long64_t entry); + + virtual void SetOption(const char *option) { fOption = option; } + virtual void SetObject(TObject *obj) { fObject = obj; } + virtual void SetInputList(TList *input) { fInput = input; } + virtual TList *GetOutputList() const { return fOutput; } + virtual void SlaveTerminate(); + virtual void Terminate(); + void CleanESD(); + void DumpSysInfo(Int_t entry); // dump system info + +protected: + TTree *fChain; //! pointer to the analyzed TTree or TChain + TTree *fTreeFriend; //! pointer to friend tree + AliESDEvent *fESDevent; //! esd event + AliESD *fESD; //! pointer to ESD + AliESDfriend *fESDfriend; //! pointer to friend + // USER defined variables + Int_t fFileNo; //! file number + TH1I *fNtracks; //! number of Tracks + TH1I *fNtracksFriend; //! number of firend Tracks + TH1I *fNClusters; //! number of clusters on track + // + // System info + // + fstream *fSysWatch; // system watch - Memory and CPU usage + fstream *fFileWatch; // file watch - write the status of the analyzed files + Int_t fDebugLevel; //debug level + + ClassDef(AliTPCSelectorESD,1); +}; + + + + + +#endif diff --git a/TPC/AliTPCSelectorTracks.cxx b/TPC/AliTPCSelectorTracks.cxx new file mode 100644 index 00000000000..8ff6a08813a --- /dev/null +++ b/TPC/AliTPCSelectorTracks.cxx @@ -0,0 +1,228 @@ +// The class definition in esdClus.h has been generated automatically +// by the ROOT utility TTree::MakeSelector(). This class is derived +// from the ROOT class TSelector. For more information on the TSelector +// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual. + +// The following methods are defined in this file: +// Begin(): called everytime a loop on the tree starts, +// a convenient place to create your histograms. +// SlaveBegin(): called after Begin(), when on PROOF called only on the +// slave servers. +// Process(): called for each event, in this function you decide what +// to read and fill your histograms. +// SlaveTerminate: called at the end of the loop on the tree, when on PROOF +// called only on the slave servers. +// Terminate(): called at the end of the loop on the tree, +// a convenient place to draw/fit your histograms. +// +// To use this file, try the following session on your Tree T: +// +// Root > T->Process("esdClus.C") +// Root > T->Process("esdClus.C","some options") +// Root > T->Process("esdClus.C+") +// +// Modification log: +// 05/11/2006 HH Correct for large pads (outer sectors) in amplitude plots + +#include "TSystem.h" +#include +#include +#include "TCint.h" +#include "TList.h" +#include "TH1I.h" +#include "TChain.h" +// +#include "AliTracker.h" +#include "AliMagF.h" +// +#include "AliESD.h" +#include "AliESDEvent.h" +#include "AliESDtrack.h" +#include "AliESDfriend.h" +#include "AliESDfriendTrack.h" +#include "AliTPCseed.h" +#include "AliTPCclusterMI.h" +// +// +#include "AliTPCcalibTracks.h" +#include "AliTPCcalibTracksGain.h" + +#include "AliTPCSelectorESD.h" +#include "AliTPCSelectorTracks.h" +#include "TProof.h" + +const char* AliTPCSelectorTracks::fgkOutputFileName = "Output.root"; + + +AliTPCSelectorTracks::AliTPCSelectorTracks(TTree *) : + AliTPCSelectorESD(), + fInit(kFALSE), + fCalibTracks(0), + fCalibTracksGain(0) +{ + // + // + // + G__SetCatchException(0); +} + +AliTPCSelectorTracks::~AliTPCSelectorTracks(){ + // + // + // +} + +void AliTPCSelectorTracks::InitComponent(){ + // + // Init Components + // + // + // USER -COMPONENT definder part + // + // before calling the process function, two objects have to be added to the chain/tree: + // clusterParam, the cluster parametrization and cuts, a AliTPCcalibTracksCuts object + // they have to be added in the following manner: + // chain->GetUserInfo()->AddLast(clusterParam); + // chain->GetUserInfo()->AddLast(cuts); + // + static Int_t counter=0; + if (!fChain){ + Error("InitComponent","ERROR - chain not initialized\n"); + } + Info("InitComponent",Form("Selector initialization No\t%d\n", counter)); + counter++; + // + + + AliTPCClusterParam *clusterParam = (AliTPCClusterParam*)fChain->GetUserInfo()->FindObject("AliTPCClusterParam"); + // + if (clusterParam == 0) Error("InitComponent","CLUSTER PARAM NOT FOUND IN CHAIN! \n"); + // + AliTPCcalibTracksCuts *cuts = (AliTPCcalibTracksCuts*)fChain->GetUserInfo()->FindObject("calibTracksCuts"); + if (cuts != 0) Info("InitComponent","cuts found in fChain! \n"); + else{ + Error("InitComponent","CUTS NOT FOUND IN CHAIN\n"); + } + if (clusterParam==0 && fInput){ + clusterParam = (AliTPCClusterParam*)fInput->FindObject("AliTPCClusterParam"); + Error("InitComponent","CLUSTER PARAM NOT FOUND IN PROOF\n"); + } + if (cuts==0 &&fInput ){ + cuts = (AliTPCcalibTracksCuts*)fInput->FindObject("calibTracksCuts"); + Error("InitComponent","CUTS NOT FOUND IN PROOF\n"); + } + if (!cuts || !clusterParam) { + if (fInput) fInput->Print(); + return; + } + + fCalibTracks = new AliTPCcalibTracks("calibTracks", "Resolution calibration object for tracks", clusterParam, cuts); + fOutput->AddLast(fCalibTracks); + + AliTPCcalibTracksGain* prevIter = (AliTPCcalibTracksGain*)fChain->GetUserInfo()->FindObject("calibTracksGain"); + if (!prevIter) Info("InitComponent", "Previous iteration of calibTracksGain not found, continuing without."); + TNamed* debugStreamPrefix = (TNamed*)fChain->GetUserInfo()->FindObject("debugStreamPrefix"); + fCalibTracksGain = new AliTPCcalibTracksGain("calibTracksGain", "Gain calibration object for tracks", cuts, debugStreamPrefix, prevIter); + fOutput->AddLast(fCalibTracksGain); + fInit=kTRUE; +} + +void AliTPCSelectorTracks::SlaveBegin(TTree * tree) +{ + // The SlaveBegin() function is called after the Begin() function. + // When running with PROOF SlaveBegin() is called on each slave server. + // The tree argument is deprecated (on PROOF 0 is passed). + + AliTPCSelectorESD::SlaveBegin(tree); + + printf(" ***** SlaveBegin ***** \n"); +} + + +void AliTPCSelectorTracks::SlaveTerminate() +{ + // The SlaveTerminate() function is called after all entries or objects + // have been processed. When running with PROOF SlaveTerminate() is called + // on each slave server. + printf ("SlaveTerminate.. \n"); + printf ("Terminate CalibTracksGain.. \n"); + if (fCalibTracksGain) fCalibTracksGain->Terminate(); +} + + + + +Int_t AliTPCSelectorTracks::ProcessIn(Long64_t entry) +{ + // + // + // + if (!fInit) InitComponent(); + if (!fInit) return 0; + Int_t status = ReadEvent(entry); + if (status<0) return status; + Int_t ntracks = (fESD) ? fESD->GetNumberOfTracks() : fESDevent->GetNumberOfTracks(); + // + // + // USER code to go here + // + AliTPCseed *seed; + + for (Int_t tr = 0; tr < ntracks; tr++){ + AliESDtrack *esdTrack = fESD ? (AliESDtrack*) fESD->GetTrack(tr): (AliESDtrack*) fESDevent->GetTrack(tr); + AliESDfriendTrack *friendtrack = (AliESDfriendTrack*) esdTrack->GetFriendTrack(); + seed = 0; + TObject *cobject = 0; + for (Int_t i = 0; ; i++){ + cobject = friendtrack->GetCalibObject(i); + if (!cobject) break; + seed = dynamic_cast(cobject); + if (seed) break; + } + + if (seed) { + fNClusters->Fill(seed->GetNumberOfClusters()); + // + fCalibTracks->Process(seed, esdTrack); // analysis is done in fCalibTracks + fCalibTracksGain->Process(seed); + } + } + CleanESD(); + return 0; +} + + +void AliTPCSelectorTracks::Terminate() +{ + // The Terminate() function is the last function to be called during + // a query. It always runs on the client, it can be used to present + // the results graphically or save the results to file. + + + if (!fOutput) return; + + TFile file(fgkOutputFileName, "recreate"); + fCalibTracksGain = (AliTPCcalibTracksGain*)fOutput->FindObject("calibTracksGain"); + // evaluate all fitters before saving them, because it doesn't seem to be possible + // to evaluate a TLinearFitter after it has been loaded from a root file + if (fCalibTracksGain) fCalibTracksGain->Evaluate(); + fOutput->Write(); + file.Close(); + printf("Successfully written file to '%s'.", fgkOutputFileName); + + + Info("Destructor","Destuctor"); + //delete fCalibTracksGain; + //delete fCalibTracks; +// printf ("Terminate... \n"); +// if (!fOutput) return; +// TFile file("Output.root","recreate"); +// printf("fOutput contains the following: \n"); +// fOutput->Print(); +// printf("Trying to write the file 'Output.root'... \n"); +// fOutput->Write(); +// file.Close(); + +} + + diff --git a/TPC/AliTPCSelectorTracks.h b/TPC/AliTPCSelectorTracks.h new file mode 100644 index 00000000000..adfb64e7fde --- /dev/null +++ b/TPC/AliTPCSelectorTracks.h @@ -0,0 +1,43 @@ + +#ifndef AliTPCSelectorTracks_h +#define AliTPCSelectorTracks_h + +#include +using namespace std; +#include + +#include +#include + + +class AliESD; +class AliESDfriend; +class TH1I; +class AliTPCcalibTracks; +class AliTPCcalibTracksGain; + + +class AliTPCSelectorTracks : public AliTPCSelectorESD { +public : + AliTPCSelectorTracks(TTree *tree=0); + virtual ~AliTPCSelectorTracks(); + virtual void SlaveBegin(TTree *tree); + virtual void SlaveTerminate(); + virtual Int_t ProcessIn(Long64_t entry); + virtual void Terminate(); + void InitComponent(); + +private: + Bool_t fInit; //! flag - component initialized + AliTPCcalibTracks *fCalibTracks; //! calib Tracks object + AliTPCcalibTracksGain *fCalibTracksGain; //! gain calibration object for tracks + static const char *fgkOutputFileName; //! filename of the output root file +// Int_t fDebugLevel; // debug level + ClassDef(AliTPCSelectorTracks,1); +}; + + + + + +#endif diff --git a/TPC/AliTPCcalibTracks.cxx b/TPC/AliTPCcalibTracks.cxx new file mode 100644 index 00000000000..05d3f9f8a20 --- /dev/null +++ b/TPC/AliTPCcalibTracks.cxx @@ -0,0 +1,2841 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + + +/////////////////////////////////////////////////////////////////////////////// +// // +// Class to analyse tracks for calibration // +// to be used as a component in AliTPCSelectorTracks // +// In the constructor you have to specify name and title // +// to get the Object out of a file. // +// The parameter 'clusterParam', a AliTPCClusterParam object // +// (needed for TPC cluster error and shape parameterization) // +// Normally you get this object out of the file 'TPCClusterParam.root' // +// In the parameter 'cuts' the cuts are specified, that decide // +// weather a track will be accepted for calibration or not. // +// // +// +// The data flow: +// +/* + Raw Data -> Local Reconstruction -> Tracking -> Calibration -> RefData (component itself) + Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam) +*/ + + // +/////////////////////////////////////////////////////////////////////////////// + +// +// ROOT includes +// +#include +#include +using namespace std; + +#include +#include +#include +#include +// +//#include +#include +#include "TLinearFitter.h" +//#include "TMatrixD.h" +#include "TTreeStream.h" +#include "TF1.h" +#include +#include +#include "TPostScript.h" +#include "TCint.h" + +#include +#include +#include +#include +#include +#include + +// +// AliROOT includes +// +#include "AliMagF.h" +#include "AliTracker.h" +#include "AliESD.h" +#include "AliESDtrack.h" +#include "AliESDfriend.h" +#include "AliESDfriendTrack.h" +#include "AliTPCseed.h" +#include "AliTPCclusterMI.h" +#include "AliTPCROC.h" + +#include "AliTPCParamSR.h" +#include "AliTrackPointArray.h" +#include "AliTPCcalibTracks.h" +#include "AliTPCClusterParam.h" +#include "AliTPCcalibTracksCuts.h" +#include "AliTPCCalPadRegion.h" +#include "AliTPCCalPad.h" +#include "AliTPCCalROC.h" +#include "TText.h" +#include "TPaveText.h" +#include "TSystem.h" + +// Thread-stuff +//#include "TThread.h" +//#include "TMutex.h" +//#include "TLockFile.h" + + +ClassImp(AliTPCcalibTracks) + + +AliTPCcalibTracks::AliTPCcalibTracks(): + TNamed() +{ + // + // AliTPCcalibTracks default constructor + // + if (fDebugLevel > 0) cout << "AliTPCcalibTracks' default constructor called" << endl; + fClusterParam = 0; +// fDebugStream = 0; + fROC = 0; + fArrayAmpRow = 0; + fArrayAmp = 0; + fArrayQDY = 0; + fArrayQDZ = 0; + fArrayQRMSY = 0; + fArrayQRMSZ = 0; + fArrayChargeVsDriftlength = 0; + fcalPadRegionChargeVsDriftlength = 0; + fDeltaY = 0; + fDeltaZ = 0; + fResolY = 0; + fResolZ = 0; + fRMSY = 0; + fRMSZ = 0; + fCuts = 0; + fHclus = 0; + fRejectedTracksHisto = 0; + fHclusterPerPadrow = 0; + fClusterCutHisto = 0; + fHclusterPerPadrowRaw = 0; + fCalPadClusterPerPadRaw = 0; + fCalPadClusterPerPad = 0; + fDebugLevel = 0; + fFitterLinY1=0; //! + fFitterLinZ1=0; //! + fFitterLinY2=0; //! + fFitterLinZ2=0; //! + fFitterParY=0; //! + fFitterParZ=0; //! + +// cout << "end of default constructor" << endl; // TO BE REMOVED +} + + +AliTPCcalibTracks::AliTPCcalibTracks(AliTPCcalibTracks* ct){ + // + // AliTPCcalibTracks copy constructor + // + if (fDebugLevel > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl; + + Bool_t dirStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + Int_t length = -1; + // backward compatibility: if the data member doesn't yet exist, it will not be merged + (ct->fArrayAmpRow) ? length = ct->fArrayAmpRow->GetEntriesFast() : length = -1; + fArrayAmpRow = new TObjArray(length); + fArrayAmp = new TObjArray(length); + for (Int_t i = 0; i < length; i++) { + fArrayAmpRow->AddAt( (TProfile*)ct->fArrayAmpRow->At(i)->Clone(), i); + fArrayAmp->AddAt( ((TProfile*)ct->fArrayAmp->At(i)->Clone()), i); + } + + (ct->fArrayQDY) ? length = ct->fArrayQDY->GetEntriesFast() : length = -1; + fArrayQDY= new TObjArray(length); + fArrayQDZ= new TObjArray(length); + fArrayQRMSY= new TObjArray(length); + fArrayQRMSZ= new TObjArray(length); + for (Int_t i = 0; i < length; i++) { + fArrayQDY->AddAt( ((TH1F*)ct->fArrayQDY->At(i)->Clone()), i); + fArrayQDZ->AddAt( ((TH1F*)ct->fArrayQDZ->At(i)->Clone()), i); + fArrayQRMSY->AddAt( ((TH1F*)ct->fArrayQRMSY->At(i)->Clone()), i); + fArrayQRMSZ->AddAt( ((TH1F*)ct->fArrayQRMSZ->At(i)->Clone()), i); + } + + (ct->fResolY) ? length = ct->fResolY->GetEntriesFast() : length = -1; + fResolY = new TObjArray(length); + fResolZ = new TObjArray(length); + fRMSY = new TObjArray(length); + fRMSZ = new TObjArray(length); + for (Int_t i = 0; i < length; i++) { + fResolY->AddAt( ((TH1F*)ct->fResolY->At(i)->Clone()), i); + fResolZ->AddAt( ((TH1F*)ct->fResolZ->At(i)->Clone()), i); + fRMSY->AddAt( ((TH1F*)ct->fRMSY->At(i)->Clone()), i); + fRMSZ->AddAt( ((TH1F*)ct->fRMSZ->At(i)->Clone()), i); + } + + (ct->fArrayChargeVsDriftlength) ? length = ct->fArrayChargeVsDriftlength->GetEntriesFast() : length = -1; + (ct->fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0; + for (Int_t i = 0; i < length; i++) { + fArrayChargeVsDriftlength->AddAt( ((TProfile*)ct->fArrayChargeVsDriftlength->At(i)->Clone()), i); + } + + fDeltaY = (TH1F*)ct->fDeltaY->Clone(); + fDeltaZ = (TH1F*)ct->fDeltaZ->Clone(); + fHclus = (TH1I*)ct->fHclus->Clone(); + fClusterCutHisto = (TH2I*)ct->fClusterCutHisto->Clone(); + fRejectedTracksHisto = (TH1I*)ct->fRejectedTracksHisto->Clone(); + fHclusterPerPadrow = (TH1I*)ct->fHclusterPerPadrow->Clone(); + fHclusterPerPadrowRaw = (TH1I*)ct->fHclusterPerPadrowRaw->Clone(); + fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)ct->fcalPadRegionChargeVsDriftlength->Clone(); + fCalPadClusterPerPad = (AliTPCCalPad*)ct->fCalPadClusterPerPad->Clone(); + fCalPadClusterPerPadRaw = (AliTPCCalPad*)ct->fCalPadClusterPerPadRaw->Clone(); + + fCuts = new AliTPCcalibTracksCuts(ct->fCuts->GetMinClusters(), ct->fCuts->GetMinRatio(), + ct->fCuts->GetMax1pt(), ct->fCuts->GetEdgeYXCutNoise(), ct->fCuts->GetEdgeThetaCutNoise()); + fDebugLevel = ct->GetLogLevel(); + SetNameTitle(ct->GetName(), ct->GetTitle()); + TH1::AddDirectory(dirStatus); // set status back to original status +// cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED +} + + +AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) : + TNamed(name, title), + fHclus(0) { + // + // AliTPCcalibTracks constructor + // specify 'name' and 'title' of your object + // specify 'clusterParam', (needed for TPC cluster error and shape parameterization) + // In the parameter 'cuts' the cuts are specified, that decide + // weather a track will be accepted for calibration or not. + // log level - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStream, 6: waste your screen + // + // All histograms are instatiated in this constructor. + // + if (fDebugLevel > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl; + G__SetCatchException(0); + + fClusterParam = clusterParam; + if (fClusterParam){ + fClusterParam->SetInstance(fClusterParam); + } + else { + 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 !!!!! + + TH1::AddDirectory(kFALSE); + + char chname[1000]; + TProfile * prof1=0; + TH1F * his1 =0; + fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160); // valgrind 3 + fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 10, 1, 10); + fHclusterPerPadrow = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160); + fHclusterPerPadrowRaw = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160); + fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad"); + fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters"); + fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159); + + // Amplitude - sector - row histograms + fArrayAmpRow = new TObjArray(72); + fArrayAmp = new TObjArray(72); + fArrayChargeVsDriftlength = new TObjArray(72); + + for (Int_t i = 0; i < 36; i++){ + sprintf(chname,"Amp_row_Sector%d",i); + prof1 = new TProfile(chname,chname,63,0,64); // valgrind 3 193,536 bytes in 354 blocks are still reachable + prof1->SetXTitle("Pad row"); + prof1->SetYTitle("Mean Max amplitude"); + fArrayAmpRow->AddAt(prof1,i); + sprintf(chname,"Amp_row_Sector%d",i+36); + prof1 = new TProfile(chname,chname,96,0,97); // valgrind 3 3,912 bytes in 6 blocks are possibly lost + prof1->SetXTitle("Pad row"); + prof1->SetYTitle("Mean Max amplitude"); + fArrayAmpRow->AddAt(prof1,i+36); + + // amplitude + sprintf(chname,"Amp_Sector%d",i); + his1 = new TH1F(chname,chname,250,0,500); // valgrind + his1->SetXTitle("Max Amplitude (ADC)"); + fArrayAmp->AddAt(his1,i); + sprintf(chname,"Amp_Sector%d",i+36); + his1 = new TH1F(chname,chname,200,0,600); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable + his1->SetXTitle("Max Amplitude (ADC)"); + fArrayAmp->AddAt(his1,i+36); + + // driftlength + sprintf(chname, "driftlengt vs. charge, ROC %i", i); + prof1 = new TProfile(chname, chname, 500, 0, 250); + prof1->SetYTitle("Charge"); + prof1->SetXTitle("Driftlength"); + fArrayChargeVsDriftlength->AddAt(prof1,i); + sprintf(chname, "driftlengt vs. charge, ROC %i", i+36); + prof1 = new TProfile(chname, chname, 500, 0, 250); + prof1->SetYTitle("Charge"); + prof1->SetXTitle("Driftlength"); + fArrayChargeVsDriftlength->AddAt(prof1,i+36); + } + + TH1::AddDirectory(kFALSE); + + fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1); + fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1); + + fResolY = new TObjArray(3); + fResolZ = new TObjArray(3); + fRMSY = new TObjArray(3); + fRMSZ = new TObjArray(3); + TH3F * his3D; + // + his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1); + fResolY->AddAt(his3D,0); + his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1); + fResolY->AddAt(his3D,1); + his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1); + fResolY->AddAt(his3D,2); + // + his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1); + fResolZ->AddAt(his3D,0); + his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1); + fResolZ->AddAt(his3D,1); + his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1); + fResolZ->AddAt(his3D,2); + // + his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8); + fRMSY->AddAt(his3D,0); + his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8); + fRMSY->AddAt(his3D,1); + his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8); + fRMSY->AddAt(his3D,2); + // + his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8); + fRMSZ->AddAt(his3D,0); + his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8); + fRMSZ->AddAt(his3D,1); + his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8); + fRMSZ->AddAt(his3D,2); + // + + TH1::AddDirectory(kFALSE); + + fArrayQDY = new TObjArray(300); + fArrayQDZ = new TObjArray(300); + fArrayQRMSY = new TObjArray(300); + fArrayQRMSZ = new TObjArray(300); + for (Int_t iq = 0; iq <= 10; iq++){ + for (Int_t ipad = 0; ipad < 3; ipad++){ + Int_t bin = GetBin(iq, ipad); + Float_t qmean = GetQ(bin); + char name[200]; + sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1); + fArrayQDY->AddAt(his3D, bin); + sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1); + fArrayQDZ->AddAt(his3D, bin); + sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1); + fArrayQRMSY->AddAt(his3D, bin); + sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean); + his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1); + fArrayQRMSZ->AddAt(his3D, bin); + } + } + + fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region"); + TProfile *tempProf; + for (UInt_t padSize = 0; padSize < 3; padSize++) { + for (UInt_t isector = 0; isector < 36; isector++) { + if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector); + if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium pads", isector); + if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long pads", isector); + tempProf = new TProfile(chname, chname, 500, 0, 250); + tempProf->SetYTitle("Charge"); + tempProf->SetXTitle("Driftlength"); + fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize); + } + } + + fFitterLinY1 = new TLinearFitter (2,"pol1"); + fFitterLinZ1 = new TLinearFitter (2,"pol1"); + fFitterLinY2 = new TLinearFitter (2,"pol1"); + fFitterLinZ2 = new TLinearFitter (2,"pol1"); + fFitterParY = new TLinearFitter (3,"pol2"); + fFitterParZ = new TLinearFitter (3,"pol2"); + + if (fDebugLevel > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; + cout << "end of main constructor" << endl; // TO BE REMOVED +} + + +AliTPCcalibTracks::~AliTPCcalibTracks() { + // + // AliTPCcalibTracks destructor + // + + if (fDebugLevel > 0) cout << "AliTPCcalibTracks' destuctor called." << endl; + Int_t length = 0; + if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast(); + for (Int_t i = 0; i < length; i++){ + delete fArrayAmpRow->At(i); + delete fArrayAmp->At(i); + } + delete fArrayAmpRow; + delete fArrayAmp; + + delete fDeltaY; + delete fDeltaZ; + + if (fResolY) length = fResolY->GetEntriesFast(); + for (Int_t i = 0; i < length; i++){ + delete fResolY->At(i); + delete fResolZ->At(i); + delete fRMSY->At(i); + delete fRMSZ->At(i); + } + delete fResolY; + delete fResolZ; + delete fRMSY; + delete fRMSZ; + + if (fArrayQDY) length = fArrayQDY->GetEntriesFast(); + for (Int_t i = 0; i < length; i++){ + delete fArrayQDY->At(i); + delete fArrayQDZ->At(i); + delete fArrayQRMSY->At(i); + delete fArrayQRMSZ->At(i); + } + + if (fArrayChargeVsDriftlength) length = fArrayChargeVsDriftlength->GetEntriesFast(); + for (Int_t i = 0; i < length; i++){ + delete fArrayChargeVsDriftlength->At(i); + } + + delete fFitterLinY1; + delete fFitterLinZ1; + delete fFitterLinY2; + delete fFitterLinZ2; + delete fFitterParY; + delete fFitterParZ; + + delete fArrayQDY; + delete fArrayQDZ; + delete fArrayQRMSY; + delete fArrayQRMSZ; + delete fArrayChargeVsDriftlength; + + delete fHclus; + delete fRejectedTracksHisto; + delete fClusterCutHisto; + delete fHclusterPerPadrow; + delete fHclusterPerPadrowRaw; + if (fCalPadClusterPerPad) delete fCalPadClusterPerPad; + if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw; + fcalPadRegionChargeVsDriftlength->Delete(); + delete fcalPadRegionChargeVsDriftlength; + if (fDebugLevel > 4) delete fDebugStream; +} + + +void AliTPCcalibTracks::AddInfo(TChain * chain, char* fileName){ + // + // Add the neccessary information for processing to the chain + // (cluster parametrization) + // + TFile clusterParamFile(fileName); + AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param"); + chain->GetUserInfo()->AddLast((TObject*)clusterParam); + cout << "Clusterparametrization added to the chain." << endl; +} + + +void AliTPCcalibTracks::Process(AliTPCseed *track, AliESDtrack *esd){ + // + // To be called in the selector + // first AcceptTrack is evaluated, then calls all the following analyse functions: + // FillResolutionHistoLocal(track) + // AlignUpDown(track, esd) + // + if (fDebugLevel > 5) Info("Process","Starting to process the track..."); + Int_t accpetStatus = AcceptTrack(track); + if (accpetStatus == 0) { + FillResolutionHistoLocal(track); + // AlignUpDown(track, esd); + } + else fRejectedTracksHisto->Fill(accpetStatus); +} + + + +Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){ + // + // calculate bins for given q and pad type + // used in TObjArray + // + Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 ); + res *= 3; + res += pad; + return res; +} + + +Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){ + // + // calculate bins for given iq and pad type + // used in TObjArray + // + return iq * 3 + pad;; +} + + +Float_t AliTPCcalibTracks::GetQ(Int_t bin){ + // + // returns to bin belonging charge + // (bin / 3 + 3)^2 + // + Int_t bin0 = bin / 3; + bin0 += 3; + return bin0 * bin0; +} + + +Float_t AliTPCcalibTracks::GetPad(Int_t bin){ + // + // returns to bin belonging pad + // bin % 3 + // + return bin % 3; +} + + + +Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){ + // + // Function, that decides wheather a given track is accepted for + // the analysis or not. + // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts' + // Returns 0 if a track is accepted or an integer different from 0 + // to indicate the failed cut + // + const Int_t kMinClusters = fCuts->GetMinClusters(); + const Float_t kMinRatio = fCuts->GetMinRatio(); + const Float_t kMax1pt = fCuts->GetMax1pt(); + const Float_t kEdgeYXCutNoise = fCuts->GetEdgeYXCutNoise(); + const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise(); + + // + // edge induced noise tracks - NEXT RELEASE will be removed during tracking + if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise ) + if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1; + if (track->GetNumberOfClusters() < kMinClusters) return 2; + Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.); + if (ratio < kMinRatio) return 3; +// Float_t mpt = track->Get1Pt(); // Get1Pt() doesn't exist any more + Float_t mpt = track->GetSigned1Pt(); + if (TMath::Abs(mpt) > kMax1pt) return 4; + //if (TMath::Abs(track->GetZ())>240.) return kFALSE; + //if (TMath::Abs(track->GetZ())<10.) return kFALSE; + //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE; + + if (fDebugLevel > 5) Info("AcceptTrack","Track has been accepted."); + return 0; +} + + +void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ + // + // fill resolution histograms - localy - tracklet in the neighborhood + // write debug information to 'TPCSelectorDebug.root' + // + // _ the main function, called during track analysis _ + // + // loop over all padrows along the track + // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction + // + // loop again over all padrows along the track + // fit tracklet (clusters in given padrow +- kDelta padrows) + // with polynom of 2nd order and two polynoms of 1st order + // take both polynoms of 1st order, calculate difference of their parameters + // add covariance matrixes and calculate chi2 of this difference + // if this chi2 is bigger than a given threshold, assume that the current cluster is + // a kink an goto next padrow + // if not: + // fill fArrayAmpRow, array with amplitudes vs. row for given sector + // fill fArrayAmp, array with amplitude histograms for give sector + // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY + // + // write debug information to 'TPCSelectorDebug.root' + // only for every kDeltaWriteDebugStream'th padrow to reduce data volume + // and to avoid redundant data + // + + if (fDebugLevel > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****"); + const Int_t kDelta = 10; // delta rows to fit + const Float_t kMinRatio = 0.75; // minimal ratio + const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal + const Float_t kErrorFraction = 0.5; // use only clusters with small interpolation error - for error param + const Int_t kFirstLargePad = 127; // medium pads -> long pads + const Float_t kLargePadSize = 1.5; // factor between medium and long pads' area + const Int_t kDeltaWriteDebugStream = 5; // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream +// TLinearFitter fFitterLinY1 = fFitterLinY1; +// TLinearFitter fFitterLinZ1 = ffFitterLinZ1; +// TLinearFitter fFitterLinY2 = ffFitterLinY2; +// TLinearFitter fFitterLinZ2 = ffFitterLinZ2; +// TLinearFitter fFitterParY = ffFitterParY; +// TLinearFitter fFitterParZ = ffFitterParZ; + TVectorD paramY0(2); + TVectorD paramZ0(2); + TVectorD paramY1(2); + TVectorD paramZ1(2); + TVectorD paramY2(3); + TVectorD paramZ2(3); + TMatrixD matrixY0(2,2); + TMatrixD matrixZ0(2,2); + TMatrixD matrixY1(2,2); + TMatrixD matrixZ1(2,2); + + // estimate mean error + Int_t nTrackletsAll = 0; // number of tracklets for given track + Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction + Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction + Int_t nClusters = 0; // working variable, number of clusters per tracklet + Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet + + fHclus->Fill(track->GetNumberOfClusters()); // for statistics overview + // --------------------------------------------------------------------- + for (Int_t irow = 0; irow < 159; irow++){ + // loop over all rows along the track + // fit tracklets (length: 13 rows) with pol2 in Y and Z direction + // calculate mean chi^2 for this track-fit in Y and Z direction + AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); + if (!cluster0) continue; // no cluster found + Int_t sector = cluster0->GetDetector(); + fHclusterPerPadrowRaw->Fill(irow); + + Int_t ipad = TMath::Nint(cluster0->GetPad()); + Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad())); + fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 ); + + if (sector != sectorG){ + // track leaves sector before it crossed enough rows to fit / initialization + nClusters = 0; + fFitterParY->ClearPoints(); + fFitterParZ->ClearPoints(); + sectorG = sector; + } + else { + nClusters++; + Double_t x = cluster0->GetX(); + fFitterParY->AddPoint(&x, cluster0->GetY(), 1); + fFitterParZ->AddPoint(&x, cluster0->GetZ(), 1); + // + if ( nClusters >= kDelta + 3 ){ + // if more than 13 (kDelta+3) clusters were added to the fitters + // fit the tracklet, increase trackletCounter + fFitterParY->Eval(); + fFitterParZ->Eval(); + nTrackletsAll++; + csigmaY += fFitterParY->GetChisquare() / (nClusters - 3.); + csigmaZ += fFitterParZ->GetChisquare() / (nClusters - 3.); + nClusters = -1; + fFitterParY->ClearPoints(); + fFitterParZ->ClearPoints(); + } + } + } // for (Int_t irow = 0; irow < 159; irow++) + // mean chi^2 for all tracklet fits in Y and in Z direction: + csigmaY = TMath::Sqrt(csigmaY / nTrackletsAll); + csigmaZ = TMath::Sqrt(csigmaZ / nTrackletsAll); + // --------------------------------------------------------------------- + + for (Int_t irow = 0; irow < 159; irow++){ + // loop again over all rows along the track + // do analysis + // + Int_t nclFound = 0; // number of clusters in the neighborhood + Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster + Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster + AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow); + if (!cluster0) continue; + Int_t sector = cluster0->GetDetector(); + Float_t xref = cluster0->GetX(); + + // Make Fit + fFitterParY->ClearPoints(); + fFitterParZ->ClearPoints(); + fFitterLinY1->ClearPoints(); + fFitterLinZ1->ClearPoints(); + fFitterLinY2->ClearPoints(); + fFitterLinZ2->ClearPoints(); + + // fit tracklet (clusters in given padrow +- kDelta padrows) + // with polynom of 2nd order and two polynoms of 1st order + // take both polynoms of 1st order, calculate difference of their parameters + // add covariance matrixes and calculate chi2 of this difference + // if this chi2 is bigger than a given threshold, assume that the current cluster is + // a kink an goto next padrow + + for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ + // loop over irow +- kDelta rows (neighboured rows) + // + // + if (idelta == 0) continue; // don't use center cluster + if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC + AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta); + if (!currentCluster) continue; + if (currentCluster->GetType() < 0) continue; + if (currentCluster->GetDetector() != sector) continue; + Double_t x = currentCluster->GetX() - xref; // x = differece: current cluster - cluster @ irow + nclFound++; + if (idelta < 0){ + ncl0++; + fFitterLinY1->AddPoint(&x, currentCluster->GetY(), csigmaY); + fFitterLinZ1->AddPoint(&x, currentCluster->GetZ(), csigmaZ); + } + if (idelta > 0){ + ncl1++; + fFitterLinY2->AddPoint(&x, currentCluster->GetY(), csigmaY); + fFitterLinZ2->AddPoint(&x, currentCluster->GetZ(), csigmaZ); + } + fFitterParY->AddPoint(&x, currentCluster->GetY(), csigmaY); + fFitterParZ->AddPoint(&x, currentCluster->GetZ(), csigmaZ); + } // loop over neighbourhood for fitter filling + + if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10); + if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow); + if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow + fFitterParY->Eval(); + fFitterParZ->Eval(); + Double_t chi2 = (fFitterParY->GetChisquare() + fFitterParZ->GetChisquare()) / (2. * nclFound - 6.); + if (chi2 > kCutChi2) fRejectedTracksHisto->Fill(9); + if (chi2 > kCutChi2) fClusterCutHisto->Fill(2, irow); + if (chi2 > kCutChi2) continue; // if chi^2 is too big goto next padrow + + // REMOVE KINK + // only when there are enough clusters (4) in each direction + if (ncl0 > 4){ + fFitterLinY1->Eval(); + fFitterLinZ1->Eval(); + } + if (ncl1 > 4){ + fFitterLinY2->Eval(); + fFitterLinZ2->Eval(); + } + + if (ncl0 > 4 && ncl1 > 4){ + fFitterLinY1->GetCovarianceMatrix(matrixY0); + fFitterLinY2->GetCovarianceMatrix(matrixY1); + fFitterLinZ1->GetCovarianceMatrix(matrixZ0); + fFitterLinZ2->GetCovarianceMatrix(matrixZ1); + fFitterLinY2->GetParameters(paramY1); + fFitterLinZ2->GetParameters(paramZ1); + fFitterLinY1->GetParameters(paramY0); + fFitterLinZ1->GetParameters(paramZ0); + paramY0 -= paramY1; + paramZ0 -= paramZ1; + matrixY0 += matrixY1; + matrixZ0 += matrixZ1; + Double_t chi2 = 0; + + TMatrixD difY(2, 1, paramY0.GetMatrixArray()); + TMatrixD difYT(1, 2, paramY0.GetMatrixArray()); + matrixY0.Invert(); + TMatrixD mulY(matrixY0, TMatrixD::kMult, difY); + TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY); + chi2 += chi2Y(0, 0); + + TMatrixD difZ(2, 1, paramZ0.GetMatrixArray()); + TMatrixD difZT(1, 2, paramZ0.GetMatrixArray()); + matrixZ0.Invert(); + TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ); + TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ); + chi2 += chi2Z(0, 0); + + // REMOVE KINK + if (chi2 * 0.25 > kCutChi2) fRejectedTracksHisto->Fill(8); + if (chi2 * 0.25 > kCutChi2) fClusterCutHisto->Fill(3, irow); + if (chi2 * 0.25 > kCutChi2) continue; // if chi2 is too big goto next padrow + // fit tracklet with polynom of 2nd order and two polynoms of 1st order + // take both polynoms of 1st order, calculate difference of their parameters + // add covariance matrixes and calculate chi2 of this difference + // if this chi2 is bigger than a given threshold, assume that the current cluster is + // a kink an goto next padrow + } + + // current padrow has no kink + + // get fit parameters from pol2 fit: + Double_t paramY[4], paramZ[4]; + paramY[0] = fFitterParY->GetParameter(0); + paramY[1] = fFitterParY->GetParameter(1); + paramY[2] = fFitterParY->GetParameter(2); + paramZ[0] = fFitterParZ->GetParameter(0); + paramZ[1] = fFitterParZ->GetParameter(1); + paramZ[2] = fFitterParZ->GetParameter(2); + + Double_t tracky = paramY[0]; + Double_t trackz = paramZ[0]; + Float_t deltay = tracky - cluster0->GetY(); + Float_t deltaz = trackz - cluster0->GetZ(); + Float_t angley = paramY[1] - paramY[0] / xref; + Float_t anglez = paramZ[1]; + + Float_t max = cluster0->GetMax(); + UInt_t isegment = cluster0->GetDetector() % 36; + Int_t padSize = 0; // short pads + if (cluster0->GetDetector() >= 36) { + padSize = 1; // medium pads + if (cluster0->GetRow() > 63) padSize = 2; // long pads + } + + // ========================================= + // wirte collected information to histograms + // ========================================= + + TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector); + if ( irow >= kFirstLargePad) max /= kLargePadSize; + profAmpRow->Fill( (Double_t)cluster0->GetRow(), max ); + TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector); + hisAmp->Fill(max); + + // remove the following two lines one day: + TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector); + profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max ); + + TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize)); + profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max ); + + fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow + Int_t ipad = TMath::Nint(cluster0->GetPad()); + Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad())); + fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 ); + + + TH3F * his3 = 0; + his3 = (TH3F*)fRMSY->At(padSize); + if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) ); + his3 = (TH3F*)fRMSZ->At(padSize); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) ); + + his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) ); + his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) ); + + + // Fill resolution histograms + Bool_t useForResol = kTRUE; + if (fFitterParY->GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE; + + if (useForResol){ + fDeltaY->Fill(deltay); + fDeltaZ->Fill(deltaz); + his3 = (TH3F*)fResolY->At(padSize); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay ); + his3 = (TH3F*)fResolZ->At(padSize); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz ); + his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay ); + his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize)); + if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz ); + } + + //============================================================================================= + + if (useForResol && nclFound > 2 * kMinRatio * kDelta + && irow % kDeltaWriteDebugStream == 0 && fDebugLevel > 4){ + if (fDebugLevel > 5) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow); + FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta); + } // if (useForResol && nclFound > 2 * kMinRatio * kDelta) + + } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++) +} // FillResolutionHistoLocal(...) + + + +void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta) { + // + // - debug part of FillResolutionHistoLocal - + // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data + // called only for fDebugLevel > 4 + // fill resolution trees + // + + Int_t sector = cluster0->GetDetector(); + Float_t xref = cluster0->GetX(); + Int_t padSize = 0; // short pads + if (cluster0->GetDetector() >= 36) { + padSize = 1; // medium pads + if (cluster0->GetRow() > 63) padSize = 2; // long pads + } + + static TLinearFitter fitY0(3, "pol2"); + static TLinearFitter fitZ0(3, "pol2"); + static TLinearFitter fitY2(5, "hyp4"); + static TLinearFitter fitZ2(5, "hyp4"); + static TLinearFitter fitY2Q(5, "hyp4"); + static TLinearFitter fitZ2Q(5, "hyp4"); + static TLinearFitter fitY2S(5, "hyp4"); + static TLinearFitter fitZ2S(5, "hyp4"); + fitY0.ClearPoints(); + fitZ0.ClearPoints(); + fitY2.ClearPoints(); + fitZ2.ClearPoints(); + fitY2Q.ClearPoints(); + fitZ2Q.ClearPoints(); + fitY2S.ClearPoints(); + fitZ2S.ClearPoints(); + + for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){ + // loop over irow +- kDelta rows (neighboured rows) + // + // + if (idelta == 0) continue; + if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC + AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta); + if (!cluster) continue; + if (cluster->GetType() < 0) continue; + if (cluster->GetDetector() != sector) continue; + Double_t x = cluster->GetX() - xref; + Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) ); + Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) ); + // + 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())), + 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 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())), + 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="<Clone(); // valgrind 3 40,139 bytes in 11 blocks are still reachable + result->SetTitle(Form("%s fit residuals",result->GetTitle())); + result->SetName(Form("%s fit residuals",result->GetName())); + TAxis *xaxis = hfit->GetXaxis(); + TAxis *yaxis = hfit->GetYaxis(); + Double_t x[2]; + for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) { + x[1] = yaxis->GetBinCenter(biny); + for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) { + x[0] = xaxis->GetBinCenter(binx); + Int_t bin = hfit->GetBin(binx, biny); + Double_t val = hfit->GetBinContent(bin); +// result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) ); + result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 ); + } + } + return result; +} + + +void AliTPCcalibTracks::SetStyle() const { + // + // set style, can be called by all draw functions + // + gROOT->SetStyle("Plain"); + gStyle->SetFillColor(10); + gStyle->SetPadColor(10); + gStyle->SetCanvasColor(10); + gStyle->SetStatColor(10); + gStyle->SetPalette(1,0); + gStyle->SetNumberContours(60); +} + + +void AliTPCcalibTracks::Draw(Option_t* opt){ + // + // draw-function of AliTPCcalibTracks + // will draws some exemplaric pictures + // + + if (fDebugLevel > 6) Info("Draw", "Drawing an exemplaric picture."); + SetStyle(); + Double_t min = 0; + Double_t max = 0; + TCanvas *c1 = new TCanvas(); + c1->Divide(0, 3); + TVirtualPad *upperThird = c1->GetPad(1); + TVirtualPad *middleThird = c1->GetPad(2); + TVirtualPad *lowerThird = c1->GetPad(3); + upperThird->Divide(2,0); + TVirtualPad *upleft = upperThird->GetPad(1); + TVirtualPad *upright = upperThird->GetPad(2); + middleThird->Divide(2,0); + TVirtualPad *middleLeft = middleThird->GetPad(1); + TVirtualPad *middleRight = middleThird->GetPad(2); + lowerThird->Divide(2,0); + TVirtualPad *downLeft = lowerThird->GetPad(1); + TVirtualPad *downRight = lowerThird->GetPad(2); + + + upleft->cd(0); + min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20; + max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20; + fDeltaY->SetAxisRange(min, max); + fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable + c1->Update(); + + upright->cd(0); + max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20; + min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20; + fDeltaZ->SetAxisRange(min, max); + fDeltaZ->Fit("gaus","q","",min, max); + c1->Update(); + + middleLeft->cd(); + fHclus->Draw(opt); + + middleRight->cd(); + fRejectedTracksHisto->Draw(opt); + TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC"); + TText *t1 = pt->AddText("1: kEdgeThetaCutNoise"); + TText *t2 = pt->AddText("2: kMinClusters"); + TText *t3 = pt->AddText("3: kMinRatio"); + TText *t4 = pt->AddText("4: kMax1pt"); + t1 = t1; t2 = t2; t3 = t3; t4 = t4; // avoid compiler warnings + pt->SetToolTipText("Legend for failed cuts"); + pt->Draw(); + + downLeft->cd(); + fHclusterPerPadrowRaw->Draw(opt); + + downRight->cd(); + fHclusterPerPadrow->Draw(opt); +} + + +void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ + // + // all functions are called, that produce pictures + // the histograms are written to the directory 'pathName' + // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file + // 'stat' is also the number of minEntries for MakeResPlotsQTree + // + + if (fDebugLevel > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName); + MakeAmpPlots(stat, pathName); + MakeDeltaPlots(pathName); + FitResolutionNew(pathName); + FitRMSNew(pathName); + MakeChargeVsDriftLengthPlots(pathName); +// MakeResPlotsQ(1, 1); + MakeResPlotsQTree(stat, pathName); +} + + +void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ + // + // creates several plots: + // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps + // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster + // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit + // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit + // Empty histograms (sectors without data) are not written to file + // the ps-files are written to the directory 'pathName', that is created if it does not exist + // 'stat': only histograms with more than 'stat' entries are written to file. + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable + TPostScript *ps; + // histograms with accumulated amplitude for all IROCs and OROCs + TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone()); + allAmpHisIROC->SetName("Amp all IROCs"); + allAmpHisIROC->SetTitle("Amp all IROCs"); + TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone()); + allAmpHisOROC->SetName("Amp all OROCs"); + allAmpHisOROC->SetTitle("Amp all OROCs"); + + ps = new TPostScript("fArrayAmp.ps", 112); + if (fDebugLevel > -1) cout << "creating fArrayAmp.ps..." << endl; + for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){ + if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat ) continue; + ps->NewPage(); + ((TH1F*)fArrayAmp->At(i))->Draw(); + c1->Update(); // valgrind 3 + if (i > 0 && i < 36) { + allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i))); + allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36))); + } + } + ps->NewPage(); + allAmpHisIROC->Draw(); + c1->Update(); // valgrind + ps->NewPage(); + allAmpHisOROC->Draw(); + c1->Update(); + ps->Close(); + delete ps; + + TH1F *his = 0; + Double_t min = 0; + Double_t max = 0; + ps = new TPostScript("fArrayAmpRow.ps", 112); + if (fDebugLevel > -1) cout << "creating fArrayAmpRow.ps..." << endl; + for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){ + his = (TH1F*)fArrayAmpRow->At(i); + if (his->GetEntries() < stat) continue; + ps->NewPage(); + min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.); + max = his->GetBinCenter(5*his->GetMaximumBin()) + 100; + his->SetAxisRange(min, max); + his->Fit("pol3", "q", "", min, max); + // his->Draw("error"); // don't use this line when you don't want to have empty pages in the ps-file + c1->Update(); + } + ps->Close(); + delete ps; + delete c1; + gSystem->ChangeDirectory(".."); +} + + +void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){ + // + // creates several plots: + // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit + // the ps-files are written to the directory 'pathName', that is created if it does not exist + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable + TPostScript *ps; + Double_t min = 0; + Double_t max = 0; + + ps = new TPostScript("DeltaYZ.ps", 112); + if (fDebugLevel > -1) cout << "creating DeltaYZ.ps..." << endl; + min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20; + max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20; + fDeltaY->SetAxisRange(min, max); + ps->NewPage(); + fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable + c1->Update(); + ps->NewPage(); + max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20; + min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20; + fDeltaZ->SetAxisRange(min, max); + fDeltaZ->Fit("gaus","q","",min, max); + c1->Update(); + ps->Close(); + delete ps; + delete c1; + gSystem->ChangeDirectory(".."); +} + + +void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){ + // + // creates charge vs. driftlength plots, one TProfile for each ROC + // is not correct like this, should be one TProfile for each sector and padsize + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable + TPostScript *ps; + ps = new TPostScript("chargeVsDriftlengthOld.ps", 112); + if (fDebugLevel > -1) cout << "creating chargeVsDriftlength.ps..." << endl; + TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone()); + TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone()); + chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC"); + chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs"); + chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC"); + chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs"); + + for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { + ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw(); + c1->Update(); + if (i > 0 && i < 36) { + chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i))); + chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36))); + } + ps->NewPage(); + } + chargeVsDriftlengthAllIROCs->Draw(); + c1->Update(); // valgrind + ps->NewPage(); + chargeVsDriftlengthAllOROCs->Draw(); + c1->Update(); + ps->Close(); + delete ps; + delete c1; + gSystem->ChangeDirectory(".."); +} + + +void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){ + // + // creates charge vs. driftlength plots, one TProfile for each ROC + // under development.... + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700)); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable +// TCanvas c1("c1", "c1", 500,(sqrt(2)*500)) + c1->Divide(0,3); + TPostScript *ps; + ps = new TPostScript("chargeVsDriftlength.ps", 111); + if (fDebugLevel > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl; + + TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone()); + TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone()); + TProfile *chargeVsDriftlengthAllLongPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone()); + chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads"); + chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads"); + chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads"); + chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads"); + chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads"); + chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads"); + + for (Int_t i = 0; i < 36; i++) { + c1->cd(1)->SetGridx(); + c1->cd(1)->SetGridy(); + ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw(); + c1->cd(2)->SetGridx(); + c1->cd(2)->SetGridy(); + ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw(); + c1->cd(3)->SetGridx(); + c1->cd(3)->SetGridy(); + ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw(); + c1->Update(); + chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)); + chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)); + chargeVsDriftlengthAllLongPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)); + ps->NewPage(); + } + c1->cd(1)->SetGridx(); + c1->cd(1)->SetGridy(); + chargeVsDriftlengthAllShortPads->Draw(); + c1->cd(2)->SetGridx(); + c1->cd(2)->SetGridy(); + chargeVsDriftlengthAllMediumPads->Draw(); + c1->cd(3)->SetGridx(); + c1->cd(3)->SetGridy(); + chargeVsDriftlengthAllLongPads->Draw(); + c1->Update(); // valgrind +// ps->NewPage(); + ps->Close(); + delete ps; + delete c1; + gSystem->ChangeDirectory(".."); +} + + + +void AliTPCcalibTracks::FitResolutionNew(char* pathName){ + // + // calculates different resulution fits in Y and Z direction + // the histograms are written to 'ResolutionYZ.ps' + // writes calculated resolution to 'resol.txt' + // all files are stored in the directory pathName + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas c; + c.Divide(2,1); + if (fDebugLevel > -1) cout << "creating ResolutionYZ.ps..." << endl; + TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112); + TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1); + fres->SetParameter(0,0.02); + fres->SetParameter(1,0.0054); + fres->SetParameter(2,0.13); + + TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory + + // create histogramw for Y-resolution + TH3F * hisResY0 = (TH3F*)fResolY->At(0); + hisResY0->FitSlicesZ(); + TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2"); + TH3F * hisResY1 = (TH3F*)fResolY->At(1); + hisResY1->FitSlicesZ(); + TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2"); + TH3F * hisResY2 = (TH3F*)fResolY->At(2); + hisResY2->FitSlicesZ(); + TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2"); + // + ps->NewPage(); + c.cd(1); + hisResY02->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost + hisResY02->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY02,fres)->Draw("surf1"); + c.Update(); + // c.SaveAs("ResolutionYPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResY12->Fit(fres, "q"); + hisResY12->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY12,fres)->Draw("surf1"); + c.Update(); + // c.SaveAs("ResolutionYPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResY22->Fit(fres, "q"); + hisResY22->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY22,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionYPad2.eps"); + + // create histogramw for Z-resolution + TH3F * hisResZ0 = (TH3F*)fResolZ->At(0); + hisResZ0->FitSlicesZ(); + TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2"); + TH3F * hisResZ1 = (TH3F*)fResolZ->At(1); + hisResZ1->FitSlicesZ(); + TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2"); + TH3F * hisResZ2 = (TH3F*)fResolZ->At(2); + hisResZ2->FitSlicesZ(); + TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2"); + + ps->NewPage(); + c.cd(1); + hisResZ02->Fit(fres, "q"); + hisResZ02->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ02,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionZPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResZ12->Fit(fres, "q"); + hisResZ12->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ12,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionZPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResZ22->Fit(fres, "q"); + hisResZ22->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ22,fres)->Draw("surf1"); + c.Update(); +// c.SaveAs("ResolutionZPad2.eps"); + ps->Close(); + delete ps; + + // write calculated resoltuions to 'resol.txt' + ofstream fresol("resol.txt"); + fresol<<"Pad 0.75 cm"<<"\n"; + hisResY02->Fit(fres, "q"); // valgrind + fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ02->Fit(fres, "q"); + fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + fresol<<"Pad 1.00 cm"<<1<<"\n"; + hisResY12->Fit(fres, "q"); // valgrind + fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ12->Fit(fres, "q"); + fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + fresol<<"Pad 1.50 cm"<<0<<"\n"; + hisResY22->Fit(fres, "q"); + fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ22->Fit(fres, "q"); + fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + + TH1::AddDirectory(kFALSE); + gSystem->ChangeDirectory(".."); + delete fres; +} + + +void AliTPCcalibTracks::FitRMSNew(char* pathName){ + // + // calculates different resulution-rms fits in Y and Z direction + // the histograms are written to 'RMS_YZ.ps' + // writes calculated resolution rms to 'rms.txt' + // all files are stored in the directory pathName + // + + SetStyle(); + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + + TCanvas c; // valgrind 3 42,120 bytes in 405 blocks are still reachable 23,816 bytes in 229 blocks are still reachable + c.Divide(2,1); + if (fDebugLevel > -1) cout << "creating RMS_YZ.ps..." << endl; + TPostScript *ps = new TPostScript("RMS_YZ.ps", 112); + TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1); + frms->SetParameter(0,0.02); + frms->SetParameter(1,0.0054); + frms->SetParameter(2,0.13); + + TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory + + // create histogramw for Y-RMS + TH3F * hisResY0 = (TH3F*)fRMSY->At(0); + hisResY0->FitSlicesZ(); + TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1"); + TH3F * hisResY1 = (TH3F*)fRMSY->At(1); + hisResY1->FitSlicesZ(); + TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1"); + TH3F * hisResY2 = (TH3F*)fRMSY->At(2); + hisResY2->FitSlicesZ(); + TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1"); + // + ps->NewPage(); + c.cd(1); + hisResY02->Fit(frms, "qn0"); + hisResY02->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY02,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSYPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResY12->Fit(frms, "qn0"); // valgrind several blocks possibly lost + hisResY12->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY12,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSYPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResY22->Fit(frms, "qn0"); + hisResY22->Draw("surf1"); + c.cd(2); + MakeDiff(hisResY22,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSYPad2.eps"); + + // create histogramw for Z-RMS + TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0); + hisResZ0->FitSlicesZ(); + TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1"); + TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1); + hisResZ1->FitSlicesZ(); + TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1"); + TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2); + hisResZ2->FitSlicesZ(); + TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1"); + // + ps->NewPage(); + c.cd(1); + hisResZ02->Fit(frms, "qn0"); // valgrind + hisResZ02->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ02,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSZPad0.eps"); + ps->NewPage(); + c.cd(1); + hisResZ12->Fit(frms, "qn0"); + hisResZ12->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ12,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSZPad1.eps"); + ps->NewPage(); + c.cd(1); + hisResZ22->Fit(frms, "qn0"); // valgrind 1 block possibly lost + hisResZ22->Draw("surf1"); + c.cd(2); + MakeDiff(hisResZ22,frms)->Draw("surf1"); + c.Update(); +// c.SaveAs("RMSZPad2.eps"); + + // write calculated resoltuion rms to 'rms.txt' + ofstream filerms("rms.txt"); + filerms<<"Pad 0.75 cm"<<"\n"; + hisResY02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost + filerms<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost + filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + filerms<<"Pad 1.00 cm"<<1<<"\n"; + hisResY12->Fit(frms, "qn0"); // valgrind 3,256 bytes in 22 blocks are indirectly lost + filerms<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ12->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable + filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + // + filerms<<"Pad 1.50 cm"<<0<<"\n"; + hisResY22->Fit(frms, "qn0"); // valgrind 40,139 bytes in 11 blocks are still reachable 330,180 bytes in 15 blocks are possibly lost + filerms<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + hisResZ22->Fit(frms, "qn0"); + filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; + + TH1::AddDirectory(kFALSE); + gSystem->ChangeDirectory(".."); + ps->Close(); + delete ps; + delete frms; +} + + +void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){ + // + // Make tree with resolution parameters + // the result is written to 'resol.root' in directory 'pathname' + // file information are available in fileInfo + // available variables in the tree 'Resol': + // Entries: number of entries for this resolution point + // nbins: number of bins that were accumulated + // Dim: direction, Dim==0: y-direction, Dim==1: z-direction + // Pad: padSize; short, medium and long + // Length: pad length, 0.75, 1, 1.5 + // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra + // Zc: center of middle bin in drift direction + // Zm: mean dirftlength for accumulated Delta-Histograms + // Zs: width of driftlength bin + // AngleC: center of middle bin in Angle-Direction + // AngleM: mean angle for accumulated Delta-Histograms + // AngleS: width of Angle-bin + // Resol: sigma for gaus fit through Delta-Histograms + // Sigma: error of sigma for gaus fit through Delta Histograms + // MeanR: mean of the Delta-Histogram + // SigmaR: rms of the Delta-Histogram + // RMSm: mean of the gaus fit through RMS-Histogram + // RMS: sigma of the gaus fit through RMS-Histogram + // RMSe0: error of mean of gaus fit in RMS-Histogram + // 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; + + gSystem->MakeDirectory(pathName); + gSystem->ChangeDirectory(pathName); + TString kFileName = "resol.root"; + 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 > -1) cout << "Histograms loaded, starting to proces..." << endl; + + //-------------------------------------------------------------------------------------------- + + char name[200]; + Double_t qMean; + Double_t zMean, angleMean, zCenter, angleCenter; + Double_t zSigma, angleSigma; + TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1); + TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1); + TF1 *fitFunction = new TF1("fitFunction", "gaus"); + Float_t entriesQ = 0; + Int_t loopCounter = 1; + + for (Int_t idim = 0; idim < 2; idim++){ + // Loop y-z corrdinate + for (Int_t ipad = 0; ipad < 3; ipad++){ + // loop pad type + for (Int_t iq = -1; iq < 10; iq++){ + // LOOP Q + if (fDebugLevel > -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; + loopCounter++; + TH3F *hres = 0; + TH3F *hrms = 0; + qMean = 0; + entriesQ = 0; + + // 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) { + hres = (TH3F*)hresl->Clone(); + hrms = (TH3F*)hrmsl->Clone(); + } + 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) { + hres = (TH3F*) hresl->Clone(); + hrms = (TH3F*) hrmsl->Clone(); + } + 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); + // 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); + // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax()); + 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; + +// projectionRes->Fit("gaus", "q0", "", xmin, xmax); +// Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2); +// Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2); + fitFunction->SetMaximum(xmax); + fitFunction->SetMinimum(xmin); + projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax); + 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; + +// projectionRms->Fit("gaus","q0","",xmin,xmax); +// Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1); +// Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2); +// Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1); +// Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2); + projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax); + 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; + + 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); + } + } // if (projectionRes->GetSum() > minEntries) + } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) + } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) + + } // iq-loop + } // ipad-loop + } // idim-loop + delete projectionRes; + delete projectionRms; + +// TFile resolFile(fTreeResol.GetFile()); + TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries)); + 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; + gSystem->ChangeDirectory(".."); +} + + + + + +/* +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 + // The object's histograms are merged via their merge functions + // 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 (!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; + collectionList->Print(); + + // create a list for each data member + TList* deltaYList = new TList; + TList* deltaZList = new TList; + TList* arrayAmpRowList = new TList; + TList* rejectedTracksList = new TList; + TList* hclusList = new TList; + TList* clusterCutHistoList = new TList; + TList* arrayAmpList = new TList; + TList* arrayQDYList = new TList; + TList* arrayQDZList = new TList; + TList* arrayQRMSYList = new TList; + TList* arrayQRMSZList = new TList; + TList* arrayChargeVsDriftlengthList = new TList; + TList* calPadRegionChargeVsDriftlengthList = new TList; + TList* hclusterPerPadrowList = new TList; + TList* hclusterPerPadrowRawList = new TList; + TList* resolYList = new TList; + TList* resolZList = new TList; + TList* rMSYList = new TList; + TList* rMSZList = new TList; + +// TList* nRowsList = new TList; +// TList* nSectList = new TList; +// TList* fileNoList = new TList; + + TIterator *listIterator = collectionList->MakeIterator(); + AliTPCcalibTracks *calibTracks = 0; + if (fDebugLevel > 1) cout << "start to iterate, filling lists" << endl; + Int_t counter = 0; + while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){ + // loop over all entries in the collectionList and get dataMembers into lists + if (!calibTracks) continue; + + deltaYList->Add( calibTracks->GetfDeltaY() ); + deltaZList->Add( calibTracks->GetfDeltaZ() ); + arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow()); + arrayAmpList->Add(calibTracks->GetfArrayAmp()); + arrayQDYList->Add(calibTracks->GetfArrayQDY()); + arrayQDZList->Add(calibTracks->GetfArrayQDZ()); + arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY()); + arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ()); + resolYList->Add(calibTracks->GetfResolY()); + resolZList->Add(calibTracks->GetfResolZ()); + rMSYList->Add(calibTracks->GetfRMSY()); + rMSZList->Add(calibTracks->GetfRMSZ()); + arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength()); + calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength()); + hclusList->Add(calibTracks->GetfHclus()); + rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto()); + clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto()); + hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow()); + hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw()); + fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad()); + fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw()); + counter++; + if (fDebugLevel > 5) cout << "filling lists, object " << counter << " added." << endl; + } + +// // reset data members +// if (fDebugLevel > 1) cout << "histogram's reset-functins are called... " << endl; +// fDeltaY->Reset(); +// fDeltaZ->Reset(); +// fHclus->Reset(); +// fClusterCutHisto->Reset(); +// fRejectedTracksHisto->Reset(); +// fHclusterPerPadrow->Reset(); +// fHclusterPerPadrowRaw->Reset(); +// for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) +// ((TProfile*)(fArrayAmpRow->At(i)))->Reset(); +// for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) +// ((TProfile*)(fArrayAmp->At(i)))->Reset(); +// for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) +// ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Reset(); +// for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) +// ((TH3F*)(fArrayQDY->At(i)))->Reset(); +// for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) +// ((TH3F*)(fArrayQDZ->At(i)))->Reset(); +// for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) +// ((TH3F*)(fArrayQRMSY->At(i)))->Reset(); +// for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) +// ((TH3F*)(fArrayQRMSZ->At(i)))->Reset(); +// for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { +// ((TH3F*)(fResolY->At(i)))->Reset(); +// ((TH3F*)(fResolZ->At(i)))->Reset(); +// ((TH3F*)(fRMSY->At(i)))->Reset(); +// ((TH3F*)(fRMSZ->At(i)))->Reset(); +// } +// for (UInt_t isec = 0; isec < 36; isec++){ +// ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(isec, 0))->Reset(); +// ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(isec, 1))->Reset(); +// ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(isec, 2))->Reset(); +// } + + // merge data members + if (fDebugLevel > 0) cout << "histogram's merge-functins are called... " << endl; + fDeltaY->Merge(deltaYList); + fDeltaZ->Merge(deltaZList); + fHclus->Merge(hclusList); + fClusterCutHisto->Merge(clusterCutHistoList); + fRejectedTracksHisto->Merge(rejectedTracksList); + fHclusterPerPadrow->Merge(hclusterPerPadrowList); + fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList); + + TObjArray* objarray = 0; + TH1* hist = 0; + TList* histList = 0; + TIterator *objListIterator = 0; + + if (fDebugLevel > 0) cout << "merging fArrayAmpRows..." << endl; + // merge fArrayAmpRows + for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) { // loop over data member, i<72 + objListIterator = arrayAmpRowList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile + if (!objarray) continue; + hist = (TProfile*)(objarray->At(i)); + histList->Add(hist); + } + ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + if (fDebugLevel > 0) cout << "merging fArrayAmps..." << endl; + // merge fArrayAmps + for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72 + TIterator *objListIterator = arrayAmpList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F + if (!objarray) continue; + hist = (TH1F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH1F*)(fArrayAmp->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + if (fDebugLevel > 0) cout << "merging fArrayQDY..." << endl; + // merge fArrayQDY + for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQDYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQDY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + if (fDebugLevel > 0) cout << "merging fArrayQDZ..." << endl; + // merge fArrayQDZ + for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQDZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + if (fDebugLevel > 0) cout << "merging fArrayQRMSY..." << endl; + // merge fArrayQRMSY + for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQRMSYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + if (fDebugLevel > 0) cout << "merging fArrayQRMSZ..." << endl; + // merge fArrayQRMSZ + for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayQRMSZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + if (fDebugLevel > 0) cout << "merging fArrayChargeVsDriftlength..." << endl; + // merge fArrayChargeVsDriftlength + for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300 + objListIterator = arrayChargeVsDriftlengthList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile + if (!objarray) continue; + hist = (TProfile*)(objarray->At(i)); + histList->Add(hist); + } + ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + if (fDebugLevel > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl; + // merge fcalPadRegionChargeVsDriftlength + AliTPCCalPadRegion *cpr = 0x0; + + /* + TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator(); + while (hist = (TProfile*)regionIterator->Next()) { + // loop over all calPadRegion's in destination calibTracks object + objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator(); + while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) { + + + hist->Merge(...); + } + */ + + for (UInt_t isec = 0; isec < 36; isec++) { + for (UInt_t padSize = 0; padSize < 3; padSize++){ + objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator(); + histList = new TList; + while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) { + // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object + if (!cpr) continue; + hist = (TProfile*)cpr->GetObject(isec, padSize); + histList->Add(hist); + } + ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList); + delete histList; + delete objListIterator; + } + } + + + + + if (fDebugLevel > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl; + // merge fResolY + for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = resolYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fResolY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + // merge fResolZ + for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = resolZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fResolZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + // merge fRMSY + for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = rMSYList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fRMSY->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + // merge fRMSZ + for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3 + objListIterator = rMSZList->MakeIterator(); + histList = new TList; + while (( objarray = (TObjArray*)objListIterator->Next() )) { + // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F + if (!objarray) continue; + hist = (TH3F*)(objarray->At(i)); + histList->Add(hist); + } + ((TH3F*)(fRMSZ->At(i)))->Merge(histList); + delete histList; + delete objListIterator; + } + + delete deltaYList; + delete deltaZList; + delete arrayAmpRowList; + delete arrayAmpList; + delete arrayQDYList; + delete arrayQDZList; + delete arrayQRMSYList; + delete arrayQRMSZList; + delete resolYList; + delete resolZList; + delete rMSYList; + delete rMSZList; +// delete nRowsList; +// delete nSectList; +// delete fileNoList; + delete listIterator; + + if (fDebugLevel > 0) cout << "merging done!" << endl; + + return 1; +} + + +AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){ + // + // function to test AliTPCcalibTrack::Merge: + // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks" + // this object is appended 'nCalTracks' times to a TList + // A new AliTPCcalibTrack object is created which merges the list + // this object is returned + // + /* + // .L AliTPCcalibTracks.cxx+g + .L libTPCcalib.so + TFile f("Output.root"); + AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks"); + //f.Close(); + TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root"); + AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param"); + clusterParamFile.Close(); + + AliTPCcalibTracks::TestMerge(calTracks, clusterParam); + */ + TList *list = new TList(); + if (ct == 0 || clusterParam == 0) return 0; + cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl; + for (Int_t i = 0; i < nCalTracks; i++) { + if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl; + list->Add(new AliTPCcalibTracks(ct)); + } + + // only for check at the end + AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(ct); + Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries(); +// Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries(); + + cout << "The list contains " << list->GetEntries() << " entries. " << endl; + + + AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018); + AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5); + cal->Merge(list); + + cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl; + cal->GetfArrayAmpRow()->At(5)->Print(); + Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries(); + + cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl; + cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " << calEntries << endl; + printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n", + calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries)); + + return cal; + +} + + + + diff --git a/TPC/AliTPCcalibTracks.h b/TPC/AliTPCcalibTracks.h new file mode 100644 index 00000000000..338fca675a5 --- /dev/null +++ b/TPC/AliTPCcalibTracks.h @@ -0,0 +1,167 @@ +#ifndef AliTPCCALIBTRACKS_H +#define AliTPCCALIBTRACKS_H + + +/////////////////////////////////////////////////////////////////////////////// +// // +// 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. // +// // +// // +/////////////////////////////////////////////////////////////////////////////// + + + +#include +class TF2; +class TH3F; +class TH1F; +class TH1I; +class TH2I; +class TH2D; +class TCollection; +class TTreeSRedirector; +class TLinearFitter; +class AliTPCClusterParam; +class TTreeSRedirector; +class AliTPCROC; +class AliTPCseed; +class AliESDtrack; +class AliTPCclusterMI; +class AliTPCcalibTracksCuts; +class AliTPCCalPadRegion; +class AliTPCCalPad; +class TChain; +class TMutex; + +using namespace std; + +class AliTPCcalibTracks : public TNamed { +public : + AliTPCcalibTracks(); // default constructor + AliTPCcalibTracks(AliTPCcalibTracks* ct); // copy constructor + AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel = 0); + virtual ~AliTPCcalibTracks(); // destructor + + static void AddInfo(TChain *chain, char *fileName); // add clusterParametrization as user info to the chain + void Process(AliTPCseed *track, AliESDtrack *esd); // to be called by the Selector + + Int_t AcceptTrack(AliTPCseed * track); + void FillResolutionHistoLocal(AliTPCseed * track); // the MAIN-FUNCTION, called for each track to fill the histograms, called by Process(...) + static TH2D* MakeDiff(TH2D * hfit, TF2 * func); + + Long64_t Merge(TCollection *li); + static AliTPCcalibTracks* TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks = 50); + + void SetStyle() const; + void Draw(Option_t* opt); // draw some exemplaric histograms for fast result-check + void MakeReport(Int_t stat, char* pathName = "plots"); // calls all functions that procude pictures, results are written to pathName, stat is the minimal statistic threshold + void MakeAmpPlots(Int_t stat, char* pathName = "plots"); + void MakeDeltaPlots(char* pathName = "plots"); + void MakeChargeVsDriftLengthPlotsOld(char* pathName = "plots"); + void MakeChargeVsDriftLengthPlots(char* pathName = "plots"); + void FitResolutionNew(char* pathName = "plots"); + void FitRMSNew(char* pathName = "plots"); + void MakeResPlotsQTree(Int_t minEntries = 100, char* pathName = "plots"); + // Thread-stuff: + // void MakeResPlotsQTreeThread(Int_t minEntries = 100, char* pathName = "plots"); + // static void* MakeResPlotsQTreeThreadFunction(void* arg); + +//protected: + TObjArray* GetfArrayAmpRow() {return fArrayAmpRow;} + TObjArray* GetfArrayAmp() {return fArrayAmp;} + TObjArray* GetfArrayQDY() {return fArrayQDY;} + TObjArray* GetfArrayQDZ() {return fArrayQDZ;} + TObjArray* GetfArrayQRMSY() {return fArrayQRMSY;} + TObjArray* GetfArrayQRMSZ() {return fArrayQRMSZ;} + TObjArray* GetfArrayChargeVsDriftlength() {return fArrayChargeVsDriftlength;} + TH1F* GetfDeltaY() {return fDeltaY;} + TH1F* GetfDeltaZ() {return fDeltaZ;} + TObjArray* GetfResolY() {return fResolY;} + TObjArray* GetfResolZ() {return fResolZ;} + TObjArray* GetfRMSY() {return fRMSY;} + TObjArray* GetfRMSZ() {return fRMSZ;} + TH1I* GetfHclus() {return fHclus;} + TH1I* GetfRejectedTracksHisto() {return fRejectedTracksHisto;} + TH1I* GetfHclusterPerPadrow() {return fHclusterPerPadrow;} + TH1I* GetfHclusterPerPadrowRaw() {return fHclusterPerPadrowRaw;} + TH2I* GetfClusterCutHisto() {return fClusterCutHisto;} + AliTPCCalPad* GetfCalPadClusterPerPad() {return fCalPadClusterPerPad; } + AliTPCCalPad* GetfCalPadClusterPerPadRaw() {return fCalPadClusterPerPadRaw;} + AliTPCCalPadRegion* GetCalPadRegionchargeVsDriftlength() {return fcalPadRegionChargeVsDriftlength;} + AliTPCcalibTracksCuts* GetCuts() {return fCuts;} + void SetLogLevel(Int_t level) {fDebugLevel = level;} + Int_t GetLogLevel() const {return fDebugLevel;} + +protected: + ClassDef(AliTPCcalibTracks,1) + + +private: + static Int_t GetBin(Float_t q, Int_t pad); + static Int_t GetBin(Int_t iq, Int_t pad); + static Float_t GetQ(Int_t bin); + static Float_t GetPad(Int_t bin); + void FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta); + AliTPCClusterParam *fClusterParam; //! pointer to cluster parameterization + TTreeSRedirector *fDebugStream; //! debug stream for + AliTPCROC *fROC; //! + TObjArray *fArrayAmpRow; // array with amplitudes versus row for given sector + TObjArray *fArrayAmp; // array with amplitude for sectors + TObjArray *fArrayQDY; // q binned delta Y histograms + TObjArray *fArrayQDZ; // q binned delta Z histograms + TObjArray *fArrayQRMSY; // q binned delta Y histograms + TObjArray *fArrayQRMSZ; // q binned delta Z histograms + TObjArray *fArrayChargeVsDriftlength; // array of arrays of TProfiles with charge vs. driftlength for each padsize and sector + AliTPCCalPadRegion *fcalPadRegionChargeVsDriftlength; // CalPadRegion, one TProfile for charge vs. driftlength for each padsize and sector + TH1F *fDeltaY; // integrated delta y histo + TH1F *fDeltaZ; // integrated delta z histo + TObjArray *fResolY; // array of resolution histograms Y + TObjArray *fResolZ; // array of resolution histograms Z + TObjArray *fRMSY; // array of RMS histograms Y + TObjArray *fRMSZ; // array of RMS histograms Z + AliTPCcalibTracksCuts *fCuts; // object with cuts, that is passed to the constructor + TH1I *fHclus; // number of clusters per track + TH1I *fRejectedTracksHisto; // histogram of rejecteced tracks, the number coresponds to the failed cut + TH1I *fHclusterPerPadrow; // histogram showing the number of clusters per padRow + TH1I *fHclusterPerPadrowRaw;// histogram showing the number of clusters per padRow before cuts on clusters are applied + TH2I *fClusterCutHisto; // histogram showing in which padRow the clusters were cutted by which criterium + AliTPCCalPad *fCalPadClusterPerPad; // AliTPCCalPad showing the number of clusters per Pad + AliTPCCalPad *fCalPadClusterPerPadRaw; // AliTPCCalPad showing the number of clusters per Pad before cuts on clusters are applied + Int_t fDebugLevel; // log level - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStream, 6: waste your screen + TLinearFitter *fFitterLinY1; //! + TLinearFitter *fFitterLinZ1; //! + TLinearFitter *fFitterLinY2; //! + TLinearFitter *fFitterLinZ2; //! + TLinearFitter *fFitterParY; //! + TLinearFitter *fFitterParZ; //! + +/* + // Threads-stuff: + static Int_t fgLoopCounter; //! only for MakeResPlotsQTreeThread, display status + static TMutex *fgWriteMutex; //! + static TMutex *fgFitResMutex; //! + static TMutex *fgFitRmsMutex; //! + struct TthreadParameterStruct { + Int_t logLevel; + Int_t minEntries; + Int_t dim; + Int_t pad; + TH3F *(*resArray)[2][3][11]; + TH3F *(*rmsArray)[2][3][11]; + TString *fileName; + TTreeSRedirector *fTreeResol; + }; + +*/ + +}; + +#endif diff --git a/TPC/AliTPCcalibTracksCuts.cxx b/TPC/AliTPCcalibTracksCuts.cxx new file mode 100644 index 00000000000..a126ce8642f --- /dev/null +++ b/TPC/AliTPCcalibTracksCuts.cxx @@ -0,0 +1,148 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + + +////////////////////////////////////////////////////// +// // +// Class to specify cuts for track analysis // +// with AliTPCcalibTracks // +// // +////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include "AliTPCcalibTracksCuts.h" + +ClassImp(AliTPCcalibTracksCuts) + +AliTPCcalibTracksCuts::AliTPCcalibTracksCuts(Int_t minClusters, Float_t minRatio, Float_t max1pt, + Float_t edgeXZCutNoise, Float_t edgeThetaCutNoise, char* outputFileName): + TNamed("calibTracksCuts", "calibTracksCuts") { + // + // Constructor for AliTPCcalibTracksCuts + // specify the cuts to be set on the processed tracks + // default cuts are for comics + // + fMinClusters = minClusters; + fMinRatio = minRatio; + fMax1pt = max1pt; + fEdgeYXCutNoise = edgeXZCutNoise; + fEdgeThetaCutNoise = edgeThetaCutNoise; + fOutputFileName = new TObjString(outputFileName); +} + +AliTPCcalibTracksCuts::~AliTPCcalibTracksCuts(){ + // + // Destructor + // + cout << "AliTPCcalibTracksCuts destructor called, nothing happend." << endl; +} + + +AliTPCcalibTracksCuts::AliTPCcalibTracksCuts(AliTPCcalibTracksCuts *cuts){ + // + // copy constructor + // + fMinClusters = cuts->GetMinClusters(); + fMinRatio = cuts->GetMinRatio(); + fMax1pt = cuts->GetMax1pt(); + fEdgeYXCutNoise = cuts->GetEdgeYXCutNoise(); + fEdgeThetaCutNoise = cuts->GetEdgeThetaCutNoise(); + fOutputFileName = new TObjString(cuts->GetOutputFileName()); +} + +AliTPCcalibTracksCuts::AliTPCcalibTracksCuts(){ + // + // default constructor + // + fMinClusters = 0; + fMinRatio = 0; + fMax1pt = 0; + fEdgeYXCutNoise = 0; + fEdgeThetaCutNoise = 0; + fOutputFileName = new TObjString(""); +} + +void AliTPCcalibTracksCuts::AddCuts(TChain * chain, char* ctype, char* outputFileName){ + // + // add predefined cuts to the chain for processing + // (creates AliTPCcalibTracksCuts object) + // the cuts are set in the following order: + // fMinClusters (number of clusters) + // fMinRatio + // fMax1pt 1 over p_t + // fEdgeYXCutNoise + // fEdgeThetaCutNoise + // + // The following predefined sets of cuts can be selected: + // laser: 20, 0.4, 0.5, 0.13, 0.018 + // cosmic: 20, 0.4, 0.5, 0.13, 0.018 + // lowflux: 20, 0.4, 5, 0.2, 0.0001 + // highflux: 20, 0.4, 5, 0.2, 0.0001 + // + + TString cutType(ctype); + cutType.ToUpper(); + AliTPCcalibTracksCuts *cuts = 0; + if (cutType == "LASER") +// cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018); + cuts = new AliTPCcalibTracksCuts(20, 0.4, 5, 0.13, 0.018, outputFileName); + else if (cutType == "COSMIC") + cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018, outputFileName); + else if (cutType == "LOWFLUX") + cuts = new AliTPCcalibTracksCuts(20, 0.4, 5, 0.2, 0.0001, outputFileName); + else if (cutType == "HIGHFLUX") + cuts = new AliTPCcalibTracksCuts(20, 0.4, 5, 0.2, 0.0001, outputFileName); + else { + cuts = new AliTPCcalibTracksCuts(20, 0.4, 5, 0.2, 0.0001, outputFileName); + cerr << "WARNING! unknown type '" << ctype << "', cuts set to default values for cosmics." << endl; + cutType = "COSMIC"; + } + chain->GetUserInfo()->AddLast(cuts); + cout << "Cuts were set to predefined set: " << cutType << endl; +} + +void AliTPCcalibTracksCuts::AddCuts(TChain * chain, Int_t minClusters, Float_t minRatio, Float_t max1pt, + Float_t edgeXZCutNoise, Float_t edgeThetaCutNoise, char* outputFileName){ + // + // add user defined cuts to the chain for processing + // (creates AliTPCcalibTracksCuts object) + // the cuts are set in the following order: + // fMinClusters (number of clusters) + // fMinRatio + // fMax1pt 1 over p_t + // fEdgeYXCutNoise + // fEdgeThetaCutNoise + // + chain->GetUserInfo()->AddLast(new AliTPCcalibTracksCuts(minClusters, minRatio, max1pt, edgeXZCutNoise, edgeThetaCutNoise, outputFileName)); + printf("Cuts were set to the individal values: minClusters: %i, minRatio: %f, max1pt: %f, edgeXZCutNoise: %f, edgeThetaCutNoise: %f \n", + minClusters, minRatio, max1pt, edgeXZCutNoise, edgeThetaCutNoise); +} + +void AliTPCcalibTracksCuts::Print(Option_t* option) { + // + // Print the cut contents + // + option = option; // to avoid compiler warnings + cout << ": The following cuts are specified: " << endl; + cout << "fMinClusters: " << fMinClusters << endl; + cout << "fMinRatio: " << fMinRatio << endl; + cout << "fMax1pt: " << fMax1pt << endl; + cout << "fEdgeYXCutNoise: " << fEdgeYXCutNoise << endl; + cout << "fEdgeThetaCutNoise: " << fEdgeThetaCutNoise << endl; + cout << "fOutputFileName: " << fOutputFileName->String().Data() << endl; +} // Prints out the specified cuts diff --git a/TPC/AliTPCcalibTracksGain.cxx b/TPC/AliTPCcalibTracksGain.cxx new file mode 100644 index 00000000000..ae158897dcf --- /dev/null +++ b/TPC/AliTPCcalibTracksGain.cxx @@ -0,0 +1,983 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//////////////////////////////////////////////////////////////////////////// +// +// === Calibration class for gain calibration using tracks === +// +// 1) Genereal idea +// ================ +// A 6-parametric parabolic function +// +// G(x, y) = p0 + p1*x + p2*y + p3*x^2 + p4*y^2 + p5 * x*y +// +// is fitted to the maximum charge values or total charge values of +// all the clusters contained in the tracks that are added to this +// object. This fit is performed for each read out chamber, in fact even +// for each type of pad sizes (thus for one segment, which consists of +// an IROC and an OROC, there are three fitters used, corresponding to +// the three pad sizes). The coordinate origin is at the center of the +// particular pad size region on each ROC. +// +// Because of the Landau nature of the charge deposition we use +// different "types" of fitters instead of one to minimize the effect +// of the long Landau tail. The difference between the fitters is only +// the charge value, that is put into them, i.e. the charge is subject +// to a transformation. At this point we use three different fit types: +// +// a) simple: the charge is put in as it is +// b) sqrt: the square root of the charge is put into the fitter +// c) log: fgkM * Log(1+q/fgkM) is put into the fitter, with +// q being the untransformed charge and fgkM=25 +// +// The results of the fits may be visualized and further used by +// creating an AliTPCCalROC or AliTPCCalPad. You may specify to undo +// the transformation and/or to normalize to the pad size. +// +// Not every track you add to this object is actually used for +// calibration. There are some cuts and conditions to exclude bad +// tracks, e.g. a pt cut to cut out tracks with too much charge +// deposition or a cut on edge clusters which are not fully +// registered and don't give a usable signal. +// +// 2) Interface / usage +// ==================== +// For each track to be added you need to call Process(). +// This method expects an AliTPCseed, which contains the necessary +// cluster information. At the moment of writing this information +// is stored in an AliESDfriend corresponding to an AliESD. +// You may also call AddTrack() if you don't want the cuts and +// other quality conditions to kick in (thus forcing the object to +// accept the track) or AddCluster() for adding single clusters. +// Call one of the Evaluate functions to evaluate the fitter(s) and +// to retrieve the fit parameters, erros and so on. You can also +// do this later on by using the different Getters. +// +// The visualization methods CreateFitCalPad() and CreateFitCalROC() +// are straight forward to use. +// +// Note: If you plan to write this object to a ROOT file, make sure +// you evaluate all the fitters *before* writing, because due +// to a bug in the fitter component writing fitters doesn't +// work properly (yet). Be aware that you cannot re-evaluate +// the fitters after loading this object from file. +// (This will be gone for a new ROOT version > v5-17-05) +// +//////////////////////////////////////////////////////////////////////////// + +#include +#include +#include "TSystem.h" +#include "TMatrixD.h" +#include "TTreeStream.h" +#include "TF1.h" +#include "AliTPCParamSR.h" +#include "AliTPCClusterParam.h" +#include "AliTrackPointArray.h" +#include "TCint.h" +#include "AliTPCcalibTracksGain.h" +#include +#include +#include +#include +#include +#include +#include + +// +// AliRoot includes +// +#include "AliMagF.h" +#include "AliMathBase.h" +// +#include "AliTPCROC.h" +#include "AliTPCParamSR.h" +#include "AliTPCCalROC.h" +#include "AliTPCCalPad.h" +// +#include "AliTracker.h" +#include "AliESD.h" +#include "AliESDtrack.h" +#include "AliESDfriend.h" +#include "AliESDfriendTrack.h" +#include "AliTPCseed.h" +#include "AliTPCclusterMI.h" +#include "AliTPCcalibTracksCuts.h" +#include "AliTPCFitPad.h" + +// REMOVE ALL OF THIS +#include +#include "AliESDEvent.h" + +/* + +TFile f("TPCCalibTracksGain.root") + +gSystem->Load("libPWG1.so") +AliTreeDraw comp +comp.SetTree(dEdx) +Double_t chi2; +TVectorD vec(3) +TMatrixD mat(3,3) +TString * str = comp.FitPlane("Cl.fQ/dedxQ.fElements[0]","Cl.fY++Cl.fX","Cl.fDetector<36",chi2,vec,mat) + +*/ + +ClassImp(AliTPCcalibTracksGain) + +const Bool_t AliTPCcalibTracksGain::fgkUseTotalCharge = kTRUE; +const Double_t AliTPCcalibTracksGain::fgkM = 25.; +const char* AliTPCcalibTracksGain::fgkDebugStreamFileName = "TPCCalibTracksGain.root"; +AliTPCParamSR* AliTPCcalibTracksGain::fgTPCparam = new AliTPCParamSR(); + +AliTPCcalibTracksGain::AliTPCcalibTracksGain() : + TNamed(), + fDebugCalPadRaw(0), + fDebugCalPadCorr(0), + fDebugStream(0), + fSimpleFitter(0), + fSqrtFitter(0), + fLogFitter(0), + fSingleSectorFitter(0), + fPrevIter(0), + fDebugStreamPrefix(0), + fCuts(0) +{ + // + // Default constructor. + // +} + +AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) : + TNamed(obj) +{ + // + // Copy constructor. + // + + fDebugCalPadRaw = new AliTPCCalPad(*(obj.fDebugCalPadRaw)); + fDebugCalPadCorr = new AliTPCCalPad(*(obj.fDebugCalPadCorr)); + fSimpleFitter = new AliTPCFitPad(*(obj.fSimpleFitter)); + fSqrtFitter = new AliTPCFitPad(*(obj.fSqrtFitter)); + fLogFitter = new AliTPCFitPad(*(obj.fLogFitter)); + fSingleSectorFitter = new AliTPCFitPad(*(obj.fSingleSectorFitter)); + fPrevIter = new AliTPCcalibTracksGain(*(obj.fPrevIter)); + fDebugStreamPrefix = new TObjString(*(obj.fDebugStreamPrefix)); + fCuts = new AliTPCcalibTracksCuts(*(obj.fCuts)); +} + +AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksGain& rhs) { + // + // Assignment operator. + // + + if (this != &rhs) { + TNamed::operator=(rhs); + fDebugCalPadRaw = new AliTPCCalPad(*(rhs.fDebugCalPadRaw)); + fDebugCalPadCorr = new AliTPCCalPad(*(rhs.fDebugCalPadCorr)); + fSimpleFitter = new AliTPCFitPad(*(rhs.fSimpleFitter)); + fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter)); + fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter)); + fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter)); + fPrevIter = new AliTPCcalibTracksGain(*(rhs.fPrevIter)); + fDebugStreamPrefix = new TObjString(*(rhs.fDebugStreamPrefix)); + fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts)); + } + return *this; +} + +AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts, TNamed* debugStreamPrefix, AliTPCcalibTracksGain* prevIter) : + TNamed(name, title), + fDebugCalPadRaw(0), + fDebugCalPadCorr(0), + fDebugStream(0), + fSimpleFitter(0), + fSqrtFitter(0), + fLogFitter(0), + fSingleSectorFitter(0), + fPrevIter(0), + fDebugStreamPrefix(0), + fCuts(0) +{ + // + // Constructor. + // + + G__SetCatchException(0); + + fCuts = cuts; + if (debugStreamPrefix) fDebugStreamPrefix = new TObjString(debugStreamPrefix->GetTitle()); + fPrevIter = prevIter; + + fSimpleFitter = new AliTPCFitPad(8, "hyp7", ""); + fSqrtFitter = new AliTPCFitPad(8, "hyp7", ""); + fLogFitter = new AliTPCFitPad(8, "hyp7", ""); + fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", ""); + + // just for debugging + fTotalTracks = 0; + fAcceptedTracks = 0; + fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction"); + fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction"); + + // this will be gone for the a new ROOT version > v5-17-05 + for (UInt_t i = 0; i < 36; i++) { + fNShortClusters[i] = 0; + fNMediumClusters[i] = 0; + fNLongClusters[i] = 0; + } + } + +AliTPCcalibTracksGain::~AliTPCcalibTracksGain() { + // + // Destructor. + // + + Info("Destructor",""); + if (fSimpleFitter) delete fSimpleFitter; + if (fSqrtFitter) delete fSqrtFitter; + if (fLogFitter) delete fLogFitter; + if (fSingleSectorFitter) delete fSingleSectorFitter; + + + if (fDebugStream) { + //fDebugStream->GetFile()->Close(); + printf("Deleting debug stream object\n"); + delete fDebugStream; + } + + if (fDebugStreamPrefix) delete fDebugStreamPrefix; + + if (fDebugCalPadRaw) delete fDebugCalPadRaw; + if (fDebugCalPadCorr) delete fDebugCalPadCorr; +} + +void AliTPCcalibTracksGain::Terminate(){ + // + // Evaluate fitters and close the debug stream. + // Also move or copy the debug stream, if a debugStreamPrefix is provided. + // + + Evaluate(); + if (fDebugStream) { + delete fDebugStream; + fDebugStream = 0; + } + + if (fDebugStreamPrefix) { + TString debugStreamPrefix = fDebugStreamPrefix->GetString(); + TString destFile(""); + destFile += debugStreamPrefix; + destFile += "/"; + destFile += gSystem->HostName(); + destFile += "_TPCCalibTracksGain.root"; + if (debugStreamPrefix.BeginsWith("root://")) { + TFile::Cp("TPCCalibTracksGain.root", destFile.Data()); + } else { + TString command("mv TPCCalibTracksGain.root "); + command += destFile; + gSystem->Exec(command.Data()); + } + } +} + +void AliTPCcalibTracksGain::AddInfo(TChain* chain, char* debugStreamPrefix, char* prevIterFileName) { + // + // Add some parameters to the chain. + // debugStreamPrefix: If specified, contains the location (either normal or xrootd directory) + // where the debug stream is moved (normal directory) or copied to (xrootd). + // prevIterFileName: If specified, contains an AliTPCcalibTracksGain object from a previous run + // for doing an iterative calibration procedure (right now unused). + // Note: The parameters are *not* added to this class, you need to do it later by retrieving + // the parameters from the chain and passing them to the constructor! + // + + if (debugStreamPrefix) { + TNamed* objDebugStreamPrefix = new TNamed("debugStreamPrefix", debugStreamPrefix); + chain->GetUserInfo()->AddLast((TObject*)objDebugStreamPrefix); + } + + if (prevIterFileName) { + TFile paramFile(prevIterFileName); + if (paramFile.IsZombie()) { + printf("File %s not found. Continuing without previous iteration.\n", prevIterFileName); + return; + } + + AliTPCcalibTracksGain *prevIter = (AliTPCcalibTracksGain*)paramFile.Get("calibTracksGain"); + if (prevIter) { + chain->GetUserInfo()->AddLast((TObject*)prevIter); + } else + printf("No calibTracksGain object found. Continuing without previous iteration.\n"); + } +} + +Bool_t AliTPCcalibTracksGain::AcceptTrack(AliTPCseed* track) { + // + // Decides whether to accept a track or not depending on track parameters and cuts + // contained as AliTPCcalibTracksCuts object fCuts. + // Tracks are discarded if the number of clusters is too low or the transverse + // momentum is too low. + // The corresponding cut values are specified in the fCuts member. + // + + if (track->GetNumberOfClusters() < fCuts->GetMinClusters()) return kFALSE; + //if ((TMath::Abs(track->GetY() / track->GetX()) > fCuts->GetEdgeYXCutNoise()) + // && (TMath::Abs(track->GetTgl()) < fCuts->GetEdgeThetaCutNoise())) return kFALSE; + //if (track->GetNumberOfClusters() / (track->GetNFoundable()+1.) < fCuts->GetMinRatio()) return kFALSE; + if (TMath::Abs(track->GetSigned1Pt()) > fCuts->GetMax1pt()) return kFALSE; + + //if (track->GetPt() < 50.) return kFALSE; + return kTRUE; +} + +void AliTPCcalibTracksGain::Process(AliTPCseed* seed) { + // + // Main method to be called when a new seed is supposed to be processed + // and be used for gain calibration. Its quality is checked before it + // is added. + // + + fTotalTracks++; + if (!AcceptTrack(seed)) return; + fAcceptedTracks++; + AddTrack(seed); +} + +Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) { + // + // Merge() merges the results of all AliTPCcalibTracksGain objects contained in + // list, thus allowing a distributed computation of several files, e.g. on PROOF. + // The merged results are merged with the data members of the AliTPCcalibTracksGain + // object used for calling the Merge method. + // The return value is 0 /*the total number of tracks used for calibration*/ if the merge + // is successful, otherwise it is -1. + // + + if (!list || list->IsEmpty()) return -1; + + if (!fSimpleFitter) fSimpleFitter = new AliTPCFitPad(8, "hyp7", ""); + if (!fSqrtFitter) fSqrtFitter = new AliTPCFitPad(8, "hyp7", ""); + if (!fLogFitter) fLogFitter = new AliTPCFitPad(8, "hyp7", ""); + if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", ""); + + // this will be gone for the a new ROOT version > v5-17-05 + for (UInt_t i = 0; i < 36; i++) { + fNShortClusters[i] = 0; + fNMediumClusters[i] = 0; + fNLongClusters[i] = 0; + } + + // just for debugging + if (fDebugCalPadRaw) delete fDebugCalPadRaw; + if (fDebugCalPadCorr) delete fDebugCalPadCorr; + fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction"); + fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction"); + fTotalTracks = 0; + fAcceptedTracks = 0; + + TIterator* iter = list->MakeIterator(); + AliTPCcalibTracksGain* cal = 0; + + while ((cal = (AliTPCcalibTracksGain*)iter->Next())) { + if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) { + Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName()); + return -1; + } + Add(cal); + } + return 0; +} + +void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) { + // + // Adds another AliTPCcalibTracksGain object to this object. + // + + fSimpleFitter->Add(cal->fSimpleFitter); + fSqrtFitter->Add(cal->fSqrtFitter); + fLogFitter->Add(cal->fLogFitter); + fSingleSectorFitter->Add(cal->fSingleSectorFitter); + + // this will be gone for the a new ROOT version > v5-17-05 + for (UInt_t iSegment = 0; iSegment < 36; iSegment++) { + fNShortClusters[iSegment] += cal->fNShortClusters[iSegment]; + fNMediumClusters[iSegment] += cal->fNMediumClusters[iSegment]; + fNLongClusters[iSegment] += cal->fNLongClusters[iSegment]; + } + + // just for debugging, remove me + fTotalTracks += cal->fTotalTracks; + fAcceptedTracks += cal->fAcceptedTracks; + fDebugCalPadRaw->Add(cal->fDebugCalPadRaw); + fDebugCalPadCorr->Add(cal->fDebugCalPadCorr); + + // Let's see later what to do with fCuts, fDebugStream and fDebugStreamPrefix, + // we'll probably won't merge them. +} + +void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) { + // + // The clusters making up the track (seed) are added to various fit functions. + // See AddCluster(...) for more detail. + // + + if (!fDebugStream) fDebugStream = new TTreeSRedirector(fgkDebugStreamFileName); + DumpTrack(seed); +} + +void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t momenta, Float_t mdedx, Int_t padType, + Float_t xcenter, TVectorD dedxQ, TVectorD dedxM, Float_t fraction, Float_t fraction2, Float_t dedge, + TVectorD parY, TVectorD parZ, TVectorD meanPos) { + // + // Adds cluster to the appropriate fitter for later analysis. + // The charge used for the fit is the maximum charge for this specific cluster or the + // accumulated charge per cluster, depending on the value of fgkUseTotalCharge. + // Depending on the pad size where the cluster is registered, the value will be put in + // the appropriate fitter. Furthermore, for each pad size three different types of fitters + // are used. The fit functions are the same for all fitters (parabolic functions), but the value + // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter + // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge + // and fgkM==25. + // + + if (!cluster) { + Error("AddCluster", "Cluster not valid."); + return; + } + + if (dedge < 3.) return; + if (fraction2 > 0.7) return; + + //Int_t padType = GetPadType(cluster->GetX()); + Double_t xx[7]; + //Double_t centerPad[2] = {0}; + //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad); + //xx[0] = cluster->GetX() - centerPad[0]; + //xx[1] = cluster->GetY() - centerPad[1]; + xx[0] = cluster->GetX() - xcenter; + xx[1] = cluster->GetY(); + xx[2] = xx[0] * xx[0]; + xx[3] = xx[1] * xx[1]; + xx[4] = xx[0] * xx[1]; + xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]); + xx[6] = xx[5] * xx[5]; + + Int_t segment = cluster->GetDetector() % 36; + Double_t q = fgkUseTotalCharge ? ((Double_t)(cluster->GetQ())) : ((Double_t)(cluster->GetMax())); // note: no normalization to pad size! + + // just for debugging + Int_t row = 0; + Int_t pad = 0; + GetRowPad(cluster->GetX(), cluster->GetY(), row, pad); + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad)); + + // correct charge by normalising to mean charge per track + q /= dedxQ[2]; + + // just for debugging + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad)); + + Double_t sqrtQ = TMath::Sqrt(q); + Double_t logQ = fgkM * TMath::Log(1 + q / fgkM); + fSimpleFitter->GetFitter(segment, padType)->AddPoint(xx, q); + fSqrtFitter->GetFitter(segment, padType)->AddPoint(xx, sqrtQ); + fLogFitter->GetFitter(segment, padType)->AddPoint(xx, logQ); + fSingleSectorFitter->GetFitter(0, padType)->AddPoint(xx, q); + + // this will be gone for the a new ROOT version > v5-17-05 + if (padType == kShortPads) + fNShortClusters[segment]++; + else if (padType == kMediumPads) + fNMediumClusters[segment]++; + else if (padType == kLongPads) + fNLongClusters[segment]++; +} + +void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) { + // + // Evaluates all fitters contained in this object. + // If the robust option is set to kTRUE a robust fit is performed with frac as + // the minimal fraction of good points (see TLinearFitter::EvalRobust for details). + // Beware: Robust fitting is much slower! + // + + fSimpleFitter->Evaluate(robust, frac); + fSqrtFitter->Evaluate(robust, frac); + fLogFitter->Evaluate(robust, frac); + fSingleSectorFitter->Evaluate(robust, frac); +} + +AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) { + // + // Creates the calibration object AliTPCcalPad using fitted parameterization + // + TObjArray tpc(72); + for (UInt_t iSector = 0; iSector < 72; iSector++) + tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize)); + return new AliTPCCalPad(&tpc); +} + +AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) { + // + // Create the AliTPCCalROC with the values per pad + // sector - sector of interest + // fitType - + // + + TVectorD par(8); + if (sector < 36) { + GetParameters(sector % 36, 0, fitType, par); + return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize); + } + else { + GetParameters(sector % 36, 1, fitType, par); + AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize); + GetParameters(sector % 36, 2, fitType, par); + AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize); + AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2); + delete roc1; + delete roc2; + return roc3; + } +} + +AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) { + // + // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the + // modifications, that the center of the region of same pad size is used as the origin + // of the fit function instead of the center of the ROC. + // The possibility of a linear fit is removed as well because it is not needed. + // Only values for pads with the given pad size are calculated, the rest is 0. + // Set undoTransformation for undoing the transformation that was applied to the + // charge values before they were put into the fitter (thus allowing comparison to the original + // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter. + // If normalizeToPadSize is true, the values are normalized to the pad size. + // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without + // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then + // applying the trafo again). + // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which + // actually doesn't describe reality! + // + + Float_t dlx, dly; + Double_t centerPad[2] = {0}; + Float_t localXY[3] = {0}; + AliTPCROC* tpcROC = AliTPCROC::Instance(); + if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector()) + return 0; + AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector); + //tpcROC->GetPositionLocal(sector, lROCfitted->GetNrows()/2, lROCfitted->GetNPads(lROCfitted->GetNrows()/2)/2, centerPad); // use this instead of the switch statement if you want to calculate the center of the ROC and not the center of the regions with the same pad size + UInt_t startRow = 0; + UInt_t endRow = 0; + switch (padType) { + case kShortPads: + startRow = 0; + endRow = lROCfitted->GetNrows(); + break; + case kMediumPads: + startRow = 0; + endRow = 64; + break; + case kLongPads: + startRow = 64; + endRow = lROCfitted->GetNrows(); + break; + } + + AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad); + Double_t value = 0; + for (UInt_t irow = startRow; irow < endRow; irow++) { + for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) { + tpcROC->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number + dlx = localXY[0] - centerPad[0]; + dly = localXY[1] - centerPad[1]; + value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly; + + // Let q' = value be the transformed value without any pad size corrections, + // let T be the transformation and let l be the pad size + // 1) don't undo transformation, don't normalize: return q' + // 2) undo transformation, don't normalize: return T^{-1} q' + // 3) undo transformation, normalize: return (T^{-1} q') / l + // 4) don't undo transformation, normalize: return T((T^{-1} q') / l) + if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */} // (1) + else { // (2), (3), (4) + //calculate T^{-1} + switch (fitType) { + case 0: /* value remains unchanged */ break; + case 1: value = value * value; break; + case 2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break; + default: Error("CreateFitCalROC", "Wrong fit type."); break; + } + if (normalizeToPadSize) value /= GetPadLength(localXY[0]); // (3) + } + if (!undoTransformation && normalizeToPadSize) { // (4) + // calculate T + switch (fitType) { + case 0: /* value remains unchanged */ break; + case 1: value = TMath::Sqrt(value); break; + case 2: value = fgkM * TMath::Log(1 + value / fgkM); break; + default: Error("CreateFitCalROC", "Wrong fit type."); break; + } + } + lROCfitted->SetValue(irow, ipad, value); + } + } + return lROCfitted; +} + +AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) { + // + // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new + // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message + // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the + // sector of the new ROC. + // + + if (!roc1 || !roc2) return 0; + if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0; + if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0; + if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch."); + AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector()); + + for (UInt_t iRow = 0; iRow < 64; iRow++) { + for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++) + roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad)); + } + for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) { + for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++) + roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad)); + } + return roc; +} + +void AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) { + // + // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType + // into the fitParam TVectorD (which should contain 8 elements). + // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter. + // Note: The fitter has to be evaluated first! + // + + GetFitter(segment, padType, fitType)->GetParameters(fitParam); +} + +void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) { + // + // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType + // into the fitError TVectorD (which should contain 8 elements). + // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter. + // Note: The fitter has to be evaluated first! + // + + GetFitter(segment, padType, fitType)->GetErrors(fitError); + fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType)); +} + +Double_t AliTPCcalibTracksGain::GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType) { + // + // Returns the reduced chi^2 value for the specified segment, padType and fitType. + // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter. + // Note: The fitter has to be evaluated first! + // + + // this will be gone for the a new ROOT version > v5-17-05 + Int_t lNClusters = 0; + switch (padType) { + case kShortPads: + lNClusters = fNShortClusters[segment]; + break; + case kMediumPads: + lNClusters = fNMediumClusters[segment]; + break; + case kLongPads: + lNClusters = fNLongClusters[segment]; + break; + } + return GetFitter(segment, padType, fitType)->GetChisquare()/(lNClusters - 8); +} + +void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) { + // + // Returns the covariance matrix for the specified segment, padType, fitType. + // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter. + // + + GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix); +} + +TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) { + // + // Returns the TLinearFitter object for the specified segment, padType, fitType. + // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter. + // + + switch (fitType) { + case kSimpleFitter: + return fSimpleFitter->GetFitter(segment, padType); + case kSqrtFitter: + return fSqrtFitter->GetFitter(segment, padType); + case kLogFitter: + return fLogFitter->GetFitter(segment, padType); + case 3: + return fSingleSectorFitter->GetFitter(0, padType); + } + return 0; +} + +Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) { + // + // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position, + // 1.5 for an OROC at long pad size position, -1 if out of bounds. + // + + Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2; + Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2; + Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2; + Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2; + Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2; + Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2; + + // if IROC + if (lx >= irocLow && lx <= irocUp) return 0.75; + // if OROC medium pads + if (lx >= orocLow1 && lx <= orocUp1) return 1.; + // if OROC long pads + if (lx >= orocLow2 && lx <= orocUp2) return 1.5; + // if out of bounds + return -1; +} + +Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) { + // + // The function returns 0 for an IROC, 1 for an OROC at medium pad size position, + // 2 for an OROC at long pad size position, -1 if out of bounds. + // + + if (GetPadLength(lx) == 0.75) return 0; + else if (GetPadLength(lx) == 1.) return 1; + else if (GetPadLength(lx) == 1.5) return 2; + return -1; +} + +// ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE +Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) { + // + // Calculate the row and pad number when the local coordinates are given. + // Returns kFALSE if the position is out of range, otherwise return kTRUE. + // WARNING: This function is preliminary and probably isn't very accurate!! + // + + Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2; + //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2; + Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2; + //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2; + Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2; + //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2; + + if (GetPadType(lx) == 0) { + row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength()); + pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth()); + } else if (GetPadType(lx) == 1) { + row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength()); + pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth()); + } else if (GetPadType(lx) == 2) { + row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength()); + pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth()); + } + else return kFALSE; + return kTRUE; +} + +void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) { + // + // Dump track information to the debug stream + // + + (*fDebugStream) << "Track" << + "Track.=" << track << // track information + "\n"; + Int_t rows[200]; + for (Int_t ipad = 0; ipad < 3; ipad++) { + GetDedx(track, ipad, rows); + } +} + +Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* rows) { + // + // GetDedx for given sector for given track + // padType - type of pads + // + + Int_t firstRow = 0, lastRow = 0; + Int_t minRow = 100; + Float_t xcenter = 0; + const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10); + const Float_t kedgey = 4.; + if (padType == 0) { + firstRow = 0; + lastRow = fgTPCparam->GetNRowLow(); + xcenter = 108.47; + } + if (padType == 1) { + firstRow = fgTPCparam->GetNRowLow(); + lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1(); + xcenter = 166.60; + } + if (padType == 2) { + firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1(); + lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp(); + xcenter = 222.6; + } + minRow = (lastRow - firstRow) / 2; + // + // + Int_t nclusters = 0; + Int_t nclustersNE = 0; // number of not edge clusters + Int_t lastSector = -1; + Float_t amplitudeQ[100]; + Float_t amplitudeM[100]; + Int_t rowIn[100]; + Int_t index[100]; + // + static TLinearFitter fitY(2, "pol1"); + static TLinearFitter fitZ(2, "pol1"); + static TVectorD parY(2); + static TVectorD parZ(2); + fitY.ClearPoints(); + fitZ.ClearPoints(); + TVectorD meanPos(6); + + for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) { + AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster); + if (cluster) { + Int_t detector = cluster->GetDetector() ; + if (lastSector == -1) lastSector = detector; + if (lastSector != detector) continue; + amplitudeQ[nclusters] = cluster->GetQ(); + amplitudeM[nclusters] = cluster->GetMax(); + rowIn[nclusters] = iCluster; + nclusters++; + Double_t dx = cluster->GetX() - xcenter; + Double_t y = cluster->GetY(); + Double_t z = cluster->GetZ(); + fitY.AddPoint(&dx, y); + fitZ.AddPoint(&dx, z); + meanPos[0] += dx; + meanPos[1] += dx; + meanPos[2] += y; + meanPos[3] += y*y; + meanPos[4] += z; + meanPos[5] += z*z; + if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++; + } + } + + if (nclusters < minRow / 2) return kFALSE; + if (nclustersNE < minRow / 2) return kFALSE; + for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters); + fitY.Eval(); + fitZ.Eval(); + fitY.GetParameters(parY); + fitZ.GetParameters(parZ); + // + // calculate truncated mean + // + TMath::Sort(nclusters, amplitudeQ, index, kFALSE); + + TVectorD dedxQ(5); + TVectorD dedxM(5); + Float_t ndedx[5]; + for (Int_t i = 0; i < 5; i++) { + dedxQ[i] = 0; + dedxM[i] = 0; + ndedx[i] = 0; + } + // + // dedx calculation + // + Int_t inonEdge = 0; + for (Int_t i = 0; i < nclusters; i++) { + Int_t rowSorted = rowIn[index[i]]; + AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted); + + if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue; //don't take edge clusters + inonEdge++; + if (inonEdge < nclustersNE * 0.5) { + ndedx[0]++; + dedxQ[0] += amplitudeQ[index[i]]; + dedxM[0] += amplitudeM[index[i]]; + } + if (inonEdge < nclustersNE * 0.6) { + ndedx[1]++; + dedxQ[1] += amplitudeQ[index[i]]; + dedxM[1] += amplitudeM[index[i]]; + } + if (inonEdge < nclustersNE * 0.7) { + ndedx[2]++; + dedxQ[2] += amplitudeQ[index[i]]; + dedxM[2] += amplitudeM[index[i]]; + } + if (inonEdge < nclustersNE * 0.8) { + ndedx[3]++; + dedxQ[3] += amplitudeQ[index[i]]; + dedxM[3] += amplitudeM[index[i]]; + } + if (inonEdge < nclustersNE * 0.9) { + ndedx[4]++; + dedxQ[4] += amplitudeQ[index[i]]; + dedxM[4] += amplitudeM[index[i]]; + } + } + for (Int_t i = 0; i < 5; i++) { + dedxQ[i] /= ndedx[i]; + dedxM[i] /= ndedx[i]; + } + + inonEdge = 0; + for (Int_t i = 0; i < nclusters; i++) { + Int_t rowSorted = rowIn[index[i]]; + AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted); + if (!cluster) { + printf("Problem\n"); + continue; + } + if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++; + Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY()); + Float_t fraction = Float_t(i) / Float_t(nclusters); + Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE); + Float_t momenta = track->GetP(); + Float_t mdedx = track->GetdEdx(); + + AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos); + + (*fDebugStream) << "dEdx" << + "Cl.=" << cluster << // cluster of interest + "P=" << momenta << // track momenta + "dedx=" << mdedx << // mean dedx - corrected for angle + "IPad=" << padType << // pad type 0..2 + "xc=" << xcenter << // x center of chamber + "dedxQ.=" << &dedxQ << // dedxQ - total charge + "dedxM.=" << &dedxM << // dedxM - maximal charge + "fraction=" << fraction << // fraction - order in statistic (0,1) + "fraction2=" << fraction2 << // fraction - order in statistic (0,1) + "dedge=" << dedge << // distance to the edge + "parY.=" << &parY << // line fit + "parZ.=" << &parZ << // line fit + "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2) + "\n"; + } + return kTRUE; +} diff --git a/TPC/AliTPCcalibTracksGain.h b/TPC/AliTPCcalibTracksGain.h new file mode 100644 index 00000000000..e80a8335b8c --- /dev/null +++ b/TPC/AliTPCcalibTracksGain.h @@ -0,0 +1,110 @@ +#ifndef AliTPCCALIBTRACKSGAIN_H +#define AliTPCCALIBTRACKSGAIN_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +#include +#include + + +#include +#include +#include +#include + +#include +#include + +using namespace std; + +class TTreeSRedirector; +class TH3F; +class TLinearFitter; + +class AliTPCClusterParam; +class AliTPCParamSR; +class AliTPCCalROC; +class AliTPCCalPad; +class AliTPCseed; +class AliTPCclusterMI; +class AliTrackPointArray; +class TTreeStream; +class AliTPCcalibTracksCuts; +class AliTPCFitPad; +class TObjString; + +class AliTPCcalibTracksGain : public TNamed { +public: + enum { + kShortPads = 0, + kMediumPads = 1, + kLongPads = 2 + }; + enum { + kSimpleFitter = 0, + kSqrtFitter = 1, + kLogFitter = 2 + }; + + AliTPCcalibTracksGain(); + AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj); + AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts, TNamed* debugStreamPrefix = 0, AliTPCcalibTracksGain* prevIter = 0); + virtual ~AliTPCcalibTracksGain(); + AliTPCcalibTracksGain& operator=(const AliTPCcalibTracksGain& rhs); + + static void AddInfo(TChain* chain, char* debugStreamPrefix = 0, char* prevIterFileName = 0); + Bool_t AcceptTrack(AliTPCseed* track); + void Terminate(); + void Add(AliTPCcalibTracksGain* cal); + void AddTrack(AliTPCseed* seed); + void AddCluster(AliTPCclusterMI* cluster, Float_t momenta, Float_t mdedx, Int_t padType, Float_t xcenter, TVectorD dedxQ, TVectorD dedxM, Float_t fraction, Float_t fraction2, Float_t dedge, TVectorD parY, TVectorD parZ, TVectorD meanPos); + void Process(AliTPCseed* seed); + Long64_t Merge(TCollection *list); + void Evaluate(Bool_t robust = kFALSE, Double_t frac = -1.); + void GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam); + void GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError); + Double_t GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType); + void GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix); + AliTPCCalPad* CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation = kFALSE, Bool_t normalizeToPadSize = kFALSE); + AliTPCCalROC* CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation = kFALSE, Bool_t normalizeToPadSize = kFALSE); + AliTPCCalROC* CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation = kFALSE, Bool_t normalizeToPadSize = kFALSE); + AliTPCCalROC* CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2); + TLinearFitter* GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType); + static Double_t GetPadLength(Double_t lx); + static Int_t GetPadType(Double_t lx); + void DumpTrack(AliTPCseed* track); + Bool_t GetDedx(AliTPCseed* track, Int_t padType, Int_t* rows); + + static Bool_t GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad); // just for debugging +protected: + + UInt_t fTotalTracks; // just for debugging + UInt_t fAcceptedTracks; // just for debugging + AliTPCCalPad* fDebugCalPadRaw; // just for debugging + AliTPCCalPad* fDebugCalPadCorr; // just for debugging + + TTreeSRedirector* fDebugStream; //! debug stream for debugging + + AliTPCFitPad* fSimpleFitter; // simple fitter for short pads + AliTPCFitPad* fSqrtFitter; // sqrt fitter for medium pads + AliTPCFitPad* fLogFitter; // log fitter for long pads + AliTPCFitPad* fSingleSectorFitter; // just for debugging + + AliTPCcalibTracksGain* fPrevIter; // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!) + + UInt_t fNShortClusters[36]; // number of clusters registered on short pads + UInt_t fNMediumClusters[36]; // number of clusters registered on medium pads + UInt_t fNLongClusters[36]; // number of clusters registered on long pads + TObjString* fDebugStreamPrefix; // pathname of the final location of the debug stream file (may also be an xrootd directory) + AliTPCcalibTracksCuts* fCuts; // cuts that are used for sieving the tracks used for calibration + + static AliTPCParamSR* fgTPCparam; //! helper object for geometry related operations + static const Double_t fgkM; // value used in the transformation of the charge values for the logarithmic fitter + static const char* fgkDebugStreamFileName; // filename of the debug stream file + static const Bool_t fgkUseTotalCharge; // whether to use the cluster's total or maximum charge + + ClassDef(AliTPCcalibTracksGain, 1); +}; + +#endif diff --git a/TPC/AliTPCcalibV0.cxx b/TPC/AliTPCcalibV0.cxx new file mode 100644 index 00000000000..3705662735a --- /dev/null +++ b/TPC/AliTPCcalibV0.cxx @@ -0,0 +1,811 @@ + +#include +#include +#include +#include +// +#include "TParticle.h" +#include "TDatabasePDG.h" +#include "AliRunLoader.h" +#include "AliStack.h" + + + +#include +#include +#include "TLinearFitter.h" +#include "TMatrixD.h" +#include "TTreeStream.h" +#include "TF1.h" + + + +#include "AliMagF.h" +#include "AliTracker.h" +#include "AliESD.h" +#include "AliESDtrack.h" +#include "AliESDfriend.h" +#include "AliESDfriendTrack.h" +#include "AliTPCseed.h" +#include "AliTPCclusterMI.h" + +#include "AliKFParticle.h" +#include "AliKFVertex.h" + +#include "AliTrackPointArray.h" +#include "TCint.h" +#include "AliTPCcalibV0.h" +#include "AliV0.h" + + + + + +ClassImp(AliTPCcalibV0) + + +AliTPCcalibV0::AliTPCcalibV0() : + TNamed(), + fDebugStream(0), + fStack(0), + fOutput(0), + fESD(0), + fPdg(0), + fParticles(0), + fV0s(0), + fGammas(0), + fV0type(0), + fTPCdEdx(0), + fTPCdEdxPi(0), + fTPCdEdxEl(0), + fTPCdEdxP(0) +{ + G__SetCatchException(0); + fDebugStream = new TTreeSRedirector("V0debug.root"); + fPdg = new TDatabasePDG; + + + // create output histograms + fTPCdEdx = new TH2F("TPCdEdX", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + BinLogX(fTPCdEdx); + fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + fTPCdEdxP = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); + + +} + +AliTPCcalibV0::~AliTPCcalibV0(){ + // + // + // + delete fDebugStream; +} + + + +void AliTPCcalibV0::ProofSlaveBegin(TList * output) +{ + // Called on PROOF - fill output list +} + + +void AliTPCcalibV0::ProcessESD(AliESD *esd, AliStack *stack){ + // + // + // + fESD = esd; + AliKFParticle::SetField(esd->GetMagneticField()); + MakeV0s(); + if (stack) { + fStack = stack; + MakeMC(); + }else{ + fStack =0; + } +} + +void AliTPCcalibV0::MakeMC(){ + // + // MC comparison + // 1. Select interesting particles + // 2. Assign the recosntructed particles + // + //1. Select interesting particles + const Float_t kMinP = 0.2; + const Float_t kMinPt = 0.1; + const Float_t kMaxR = 0.5; + const Float_t kMaxTan = 1.2; + const Float_t kMaxRad = 150; + // + if (!fParticles) fParticles = new TObjArray; + TParticle *part=0; + // + Int_t entries = fStack->GetNtrack(); + for (Int_t ipart=0; ipartParticle(ipart); + if (!part) continue; + if (part->P()R()>kMaxR) continue; + if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue; + Bool_t isInteresting = kFALSE; + if (part->GetPdgCode()==22) isInteresting =kTRUE; + if (part->GetPdgCode()==310) isInteresting =kTRUE; + if (part->GetPdgCode()==111) isInteresting =kTRUE; + if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE; + + // + if (!isInteresting) continue; + fParticles->AddLast(new TParticle(*part)); + } + if (fParticles->GetEntries()<1) { + return; + } + // + // + // + Int_t sentries=fParticles->GetEntries();; + for (Int_t ipart=0; ipartAt(ipart); + TParticle *p0 = 0; + TParticle *p1 = 0; + + Int_t nold =0; + Int_t nnew =0; + Int_t id0 = part->GetDaughter(0); + Int_t id1 = part->GetDaughter(1); + if (id0>=fStack->GetNtrack() ) id0*=-1; + if (id1>=fStack->GetNtrack() ) id1*=-1; + Bool_t findable = kTRUE; + if (id0<0 || id1<0) findable = kFALSE; + Int_t charge =0; + if (findable){ + p0 = fStack->Particle(id0); + if (p0->R()>kMaxRad) findable = kFALSE; + if (p0->Pt()Vz()>250) findable= kFALSE; + if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; + if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE; + else + if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++; + + p1 = fStack->Particle(id1); + if (p1->R()>kMaxRad) findable = kFALSE; + if (p1->Pt()Vz())>250) findable= kFALSE; + if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; + if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE; + else + if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++; + + } + // (*fDebugStream)<<"MC0"<< + // "P.="<Pt(), p1->Pt()); + Int_t type=-1; + + // + // + AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); + for (Int_t ivertex=0; ivertexGetNumberOfV0s(); ivertex++){ + AliESDv0 * v0 = fESD->GetV0(ivertex); + // select coresponding track + AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); + if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue; + AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); + if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue; + TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel())); + TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel())); + // + // + if ( v0->GetOnFlyStatus()) nnew++; + if (!v0->GetOnFlyStatus()) nold++; + if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 && TMath::Abs(pp->GetPdgCode())==11) + type =1; + if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 && TMath::Abs(pp->GetPdgCode())==211) + type =0; + if (part->GetPdgCode()==3122){ + if (TMath::Abs(pn->GetPdgCode())==210 ) type=2; + else type=3; + } + AliKFParticle *v0kf = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); + v0kf->SetProductionVertex( primVtx ); + AliKFParticle *v0kfc = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); + v0kfc->SetProductionVertex( primVtx ); + v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass()); + Float_t chi2 = v0kf->GetChi2(); + Float_t chi2C = v0kf->GetChi2(); + // + // + (*fDebugStream)<<"MCRC"<< + "P.="<Delete(); + +} + + + +void AliTPCcalibV0::MakeV0s(){ + // + // + // + const Int_t kMinCluster=40; + const Float_t kMinR =0; + if (! fV0s) fV0s = new TObjArray(10); + fV0s->Clear(); + // + // Old V0 finder + // + for (Int_t ivertex=0; ivertexGetNumberOfV0s(); ivertex++){ + AliESDv0 * v0 = fESD->GetV0(ivertex); + if (v0->GetOnFlyStatus()) continue; + fV0s->AddLast(v0); + } + ProcessV0(0); + fV0s->Clear(0); + // + // MI V0 finder + // + for (Int_t ivertex=0; ivertexGetNumberOfV0s(); ivertex++){ + AliESDv0 * v0 = fESD->GetV0(ivertex); + if (!v0->GetOnFlyStatus()) continue; + fV0s->AddLast(v0); + } + ProcessV0(1); + fV0s->Clear(); + return; + // + // combinatorial + // + Int_t ntracks = fESD->GetNumberOfTracks(); + for (Int_t itrack0=0; itrack0GetTrack(itrack0); + if (track0->GetSign()>0) continue; + if ( track0->GetTPCNcls()GetKinkIndex(0)>0) continue; + // + for (Int_t itrack1=0; itrack1GetTrack(itrack1); + if (track1->GetSign()<0) continue; + if ( track1->GetTPCNcls()GetKinkIndex(0)>0) continue; + // + // AliExternalTrackParam param0(*track0); + // AliExternalTrackParam param1(*track1); + AliV0 vertex; + vertex.SetParamN(*track0); + vertex.SetParamP(*track1); + Float_t xyz[3]; + xyz[0] = fESD->GetPrimaryVertex()->GetXv(); + xyz[1] = fESD->GetPrimaryVertex()->GetYv(); + xyz[2] = fESD->GetPrimaryVertex()->GetZv(); + vertex.Update(xyz); + if (vertex.GetRr()1.) continue; + if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue; + // if (vertex.GetPointAngle()<0.9) continue; + vertex.SetIndex(0,itrack0); + vertex.SetIndex(1,itrack1); + fV0s->AddLast(new AliV0(vertex)); + } + } + ProcessV0(2); + for (Int_t i=0;iGetEntries(); i++) delete fV0s->At(i); + fV0s->Clear(); +} + + + + + + +// void AliTPCcalibV0::ProcessV0(Int_t ftype){ +// // +// // +// const Double_t ktimeK0 = 2.684; +// const Double_t ktimeLambda = 7.89; + + +// if (! fGammas) fGammas = new TObjArray(10); +// fGammas->Clear(); +// Int_t nV0s = fV0s->GetEntries(); +// if (nV0s==0) return; +// AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); +// // +// for (Int_t ivertex=0; ivertexAt(ivertex); +// AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); +// AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); +// // +// // +// // +// AliKFParticle *v0K0 = Fit(primVtx,v0,211,211); +// AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11); +// AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211); +// AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212); +// //Set production vertex +// v0K0->SetProductionVertex( primVtx ); +// v0Gamma->SetProductionVertex( primVtx ); +// v0Lambda42->SetProductionVertex( primVtx ); +// v0Lambda24->SetProductionVertex( primVtx ); +// Double_t massK0, massGamma, massLambda42,massLambda24, massSigma; +// v0K0->GetMass( massK0,massSigma); +// v0Gamma->GetMass( massGamma,massSigma); +// v0Lambda42->GetMass( massLambda42,massSigma); +// v0Lambda24->GetMass( massLambda24,massSigma); +// Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF(); +// Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF(); +// Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF(); +// Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF(); +// // +// // Mass Contrained params +// // +// AliKFParticle *v0K0C = Fit(primVtx,v0,211,211); +// AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); +// AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211); +// AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212); +// // +// v0K0C->SetProductionVertex( primVtx ); +// v0GammaC->SetProductionVertex( primVtx ); +// v0Lambda42C->SetProductionVertex( primVtx ); +// v0Lambda24C->SetProductionVertex( primVtx ); + +// v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); +// v0GammaC->SetMassConstraint(0); +// v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); +// v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); +// // +// Double_t timeK0, sigmaTimeK0; +// Double_t timeLambda42, sigmaTimeLambda42; +// Double_t timeLambda24, sigmaTimeLambda24; +// v0K0C->GetLifeTime(timeK0, sigmaTimeK0); +// //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0); +// v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42); +// v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24); + + +// // +// Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF(); +// if (chi2K0C<0) chi2K0C=100; +// Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF(); +// if (chi2GammaC<0) chi2GammaC=100; +// Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF(); +// if (chi2Lambda42C<0) chi2Lambda42C=100; +// Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF(); +// if (chi2Lambda24C<0) chi2Lambda24C=100; +// // +// Float_t minChi2C=99; +// Int_t type =-1; +// if (chi2K0CGetP()/fPdg->GetParticle(-2212)->Mass(); +// Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass(); +// Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass(); +// Float_t dedxTeorP = TPCBetheBloch(betaGammaP); +// Float_t dedxTeorPi = TPCBetheBloch(betaGammaPi);; +// Float_t dedxTeorEl = TPCBetheBloch(betaGammaEl);; +// // +// // +// if (minChi2>50) continue; +// (*fDebugStream)<<"V0"<< +// "ftype="<AddLast(v0); +// // +// // +// // +// delete v0K0; +// delete v0Gamma; +// delete v0Lambda42; +// delete v0Lambda24; +// delete v0K0C; +// delete v0GammaC; +// delete v0Lambda42C; +// delete v0Lambda24C; +// } +// ProcessPI0(); +// } + + + + +void AliTPCcalibV0::ProcessV0(Int_t ftype){ + // + // + + if (! fGammas) fGammas = new TObjArray(10); + fGammas->Clear(); + Int_t nV0s = fV0s->GetEntries(); + if (nV0s==0) return; + AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); + // + for (Int_t ivertex=0; ivertexAt(ivertex); + AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track + AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track + + const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam(); + const AliExternalTrackParam * paramInPos = trackP->GetInnerParam(); + + if (!paramInPos) continue; // in case the inner paramters do not exist + if (!paramInNeg) continue; + // + // + // + AliKFParticle *v0K0 = Fit(primVtx,v0,-211,211); + AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11); + AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211); + AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212); + //Set production vertex + v0K0->SetProductionVertex( primVtx ); + v0Gamma->SetProductionVertex( primVtx ); + v0Lambda42->SetProductionVertex( primVtx ); + v0Lambda24->SetProductionVertex( primVtx ); + Double_t massK0, massGamma, massLambda42,massLambda24, massSigma; + v0K0->GetMass( massK0,massSigma); + v0Gamma->GetMass( massGamma,massSigma); + v0Lambda42->GetMass( massLambda42,massSigma); + v0Lambda24->GetMass( massLambda24,massSigma); + Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF(); + Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF(); + Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF(); + Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF(); + // + // Mass Contrained params + // + AliKFParticle *v0K0C = Fit(primVtx,v0,-211,211); + AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); + AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar + AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda + // + v0K0C->SetProductionVertex( primVtx ); + v0GammaC->SetProductionVertex( primVtx ); + v0Lambda42C->SetProductionVertex( primVtx ); + v0Lambda24C->SetProductionVertex( primVtx ); + + v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); + v0GammaC->SetMassConstraint(0); + v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); + v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); + // + Double_t timeK0, sigmaTimeK0; + Double_t timeLambda42, sigmaTimeLambda42; + Double_t timeLambda24, sigmaTimeLambda24; + v0K0C->GetLifeTime(timeK0, sigmaTimeK0); + //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0); + v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42); + v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24); + + + // + Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF(); + if (chi2K0C<0) chi2K0C=100; + Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF(); + if (chi2GammaC<0) chi2GammaC=100; + Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF(); + if (chi2Lambda42C<0) chi2Lambda42C=100; + Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF(); + if (chi2Lambda24C<0) chi2Lambda24C=100; + // + Float_t minChi2C=99; + Int_t type =-1; + if (chi2K0CGetP()/fPdg->GetParticle(-211)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); + break; + case 1: + betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass(); + break; + case 2: + betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); + break; + case 3: + betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); + betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass(); + break; + } + + // cuts and histogram filling + Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2 + + if (minChi2C < 2 && ftype == 1) { + // + if (chi2K0C < 10*minChi2C) numCand++; + if (chi2GammaC < 10*minChi2C) numCand++; + if (chi2Lambda42C < 10*minChi2C) numCand++; + if (chi2Lambda24C < 10*minChi2C) numCand++; + // + if (numCand < 2) { + if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal()); + if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal()); + } + } + + // + // + // write output tree + if (minChi2>50) continue; + (*fDebugStream)<<"V0"<< + "ftype="<AddLast(v0); + // + // + // + delete v0K0; + delete v0Gamma; + delete v0Lambda42; + delete v0Lambda24; + delete v0K0C; + delete v0GammaC; + delete v0Lambda42C; + delete v0Lambda24C; + } + ProcessPI0(); +} + + + +void AliTPCcalibV0::ProcessPI0(){ + // + // + // + Int_t nentries = fGammas->GetEntries(); + if (nentries<2) return; + // + Double_t m0[3], m1[3]; + AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); + for (Int_t i0=0; i0At(i0); + v00->GetPxPyPz (m0[0], m0[1], m0[2]); + AliKFParticle *p00 = Fit(primVtx, v00, 11,-11); + p00->SetProductionVertex( primVtx ); + p00->SetMassConstraint(0); + // + for (Int_t i1=i0; i1At(i1); + v01->GetPxPyPz (m1[0], m1[1], m1[2]); + AliKFParticle *p01 = Fit(primVtx, v01, 11,-11); + p01->SetProductionVertex( primVtx ); + p01->SetMassConstraint(0); + if (v00->GetIndex(0) != v01->GetIndex(0) && + v00->GetIndex(1) != v01->GetIndex(1)){ + AliKFParticle pi0( *p00,*p01); + pi0.SetProductionVertex(primVtx); + Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]); + Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]); + Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2]))); + (*fDebugStream)<<"PI0"<< + "v00.="<GetParamN()), PDG1 ); + AliKFParticle p2( *(v0->GetParamP()), PDG2 ); + AliKFParticle *V0 = new AliKFParticle; + Double_t x, y, z; + v0->GetXYZ(x,y,z ); + V0->SetVtxGuess(x,y,z); + *(V0)+=p1; + *(V0)+=p2; + return V0; +} + + + +Float_t AliTPCcalibV0::TPCBetheBloch(Float_t bg) +{ + // + // Bethe-Bloch energy loss formula + // + const Double_t kp1=0.76176e-1; + const Double_t kp2=10.632; + const Double_t kp3=0.13279e-4; + const Double_t kp4=1.8631; + const Double_t kp5=1.9479; + Double_t dbg = (Double_t) bg; + Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg); + Double_t aa = TMath::Power(beta,kp4); + Double_t bb = TMath::Power(1./dbg,kp5); + bb=TMath::Log(kp3+bb); + return ((Float_t)((kp2-aa-bb)*kp1/aa)); +} + + + + + + + +TH2F * AliTPCcalibV0::GetHistograms() { + // + // + // + return fTPCdEdx; +} + + + + +void AliTPCcalibV0::BinLogX(TH2F *h) { + // + // + // + TAxis *axis = h->GetXaxis(); + int bins = axis->GetNbins(); + + Double_t from = axis->GetXmin(); + Double_t to = axis->GetXmax(); + Double_t *new_bins = new Double_t[bins + 1]; + + new_bins[0] = from; + Double_t factor = pow(to/from, 1./bins); + + for (int i = 1; i <= bins; i++) { + new_bins[i] = factor * new_bins[i-1]; + } + axis->Set(bins, new_bins); + delete new_bins; + +} + + + + diff --git a/TPC/AliTPCcalibV0.h b/TPC/AliTPCcalibV0.h new file mode 100644 index 00000000000..8004f76de1c --- /dev/null +++ b/TPC/AliTPCcalibV0.h @@ -0,0 +1,66 @@ +#ifndef AliTPCCALIBV0_H +#define AliTPCCALIBV0_H + + +#include + + +class TTreeSRedirector; +class AliTPCROC; +class AliTPCseed; +class AliESDtrack; +class AliESD; +class TH3F; +class TH1F; +class TH2F; +class TH1I; +class TDatabasePDG; +class AliKFParticle; +class AliKFVertex; +class AliESDv0; +class TArrayI; +class TTree; +class AliStack; + +class AliTPCcalibV0 : public TNamed { +public : + + // List of branches + + AliTPCcalibV0(); + virtual ~AliTPCcalibV0(); + virtual void ProofSlaveBegin(TList * output); + void ProcessESD(AliESD *esd, AliStack *stack=0); + void MakeMC(); + void MakeV0s(); + void ProcessV0(Int_t ftype); + void ProcessPI0(); + Float_t TPCBetheBloch(Float_t bg); + TH2F * GetHistograms(); + void BinLogX(TH2F * h); + // + // + // +protected: + AliKFParticle * Fit(AliKFVertex & primVtx, AliESDv0 *v0, Int_t PDG1, Int_t PDG2); +private: + TTreeSRedirector *fDebugStream; //debug stream for + AliStack *fStack; // pointer to kinematic tree + TList *fOutput; //output list + AliESD *fESD; //! current ED to proccess - NOT OWNER + TDatabasePDG *fPdg; // particle database + TObjArray *fParticles; // array of selected MC particles + TObjArray *fV0s; // array of V0s + TObjArray *fGammas; // gamma conversion candidates + // + TArrayI *fV0type; // array of types for V0s + TH2F *fTPCdEdx; // dEdx spectra + TH2F *fTPCdEdxPi; // dEdx spectra - pion anti-pion + TH2F *fTPCdEdxEl; // dEdx spectra - electroen -positrons from gamma + TH2F *fTPCdEdxP; // dEdx spectra - proton antiproton - lambda - antilambda + // + ClassDef(AliTPCcalibV0,1); +}; + + +#endif diff --git a/TPC/TPCcalibLinkDef.h b/TPC/TPCcalibLinkDef.h new file mode 100644 index 00000000000..f78fcfe5142 --- /dev/null +++ b/TPC/TPCcalibLinkDef.h @@ -0,0 +1,28 @@ +#ifdef __CINT__ +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + + + +#pragma link C++ class AliTPCcalibTracks+; +#pragma link C++ class AliTPCcalibTracksCuts+; +#pragma link C++ class AliTPCcalibTracksGain+; +#pragma link C++ class AliTPCFitPad+; +#pragma link C++ class AliTPCCalPadRegion+; + +#pragma link C++ class AliTPCSelectorESD+; +#pragma link C++ class AliTPCSelectorTracks+; + + +#endif + + + + + diff --git a/TPC/libTPCcalib.pkg b/TPC/libTPCcalib.pkg new file mode 100644 index 00000000000..ac2d42b20f2 --- /dev/null +++ b/TPC/libTPCcalib.pkg @@ -0,0 +1,14 @@ + +SRCS = AliTPCcalibTracksCuts.cxx AliTPCcalibTracks.cxx AliTPCcalibTracksGain.cxx \ + AliTPCSelectorESD.cxx AliTPCSelectorTracks.cxx AliTPCCalPadRegion.cxx AliTPCFitPad.cxx + + +HDRS:= $(SRCS:.cxx=.h) + +EINCLUDE:=RAW TPC STEER + +DHDR= TPCcalibLinkDef.h + +#EXPORT:=$(SRCS:.cxx=.h) + + -- 2.43.5