#include <iostream>
+#include <fstream>
#include <TDatabasePDG.h>
#include <TRandom.h>
#include <TSpline.h>
#include <TObjString.h>
#include <TROOT.h>
+#include <TSystem.h>
+#include <TObjArray.h>
#include <AliLog.h>
#include <AliTPCROC.h>
,fEvent(0x0)
,fCurrentTrack(0)
,fTPCCorrection(0x0)
+ ,fSCList(0x0)
+ ,fSCListFile()
,fCorrectionFile("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.lookup.root")
,fOutputFileName("toyMC.root")
,fOutFile(0x0)
,fUseStepCorrection(kFALSE)
,fUseMaterialBudget(kFALSE)
,fIsLaser(kTRUE)
+ ,fPrereadSCList(kFALSE)
{
fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
fTPCParam->ReadGeoMatrices();
,fEvent(0x0)
,fCurrentTrack(0)
,fTPCCorrection(gen.fTPCCorrection)
+ ,fSCList(gen.fSCList)
+ ,fSCListFile(gen.fSCListFile)
,fCorrectionFile(gen.fCorrectionFile)
,fOutputFileName(gen.fOutputFileName)
,fOutFile(0x0)
,fUseStepCorrection(gen.fUseStepCorrection)
,fUseMaterialBudget(gen.fUseMaterialBudget)
,fIsLaser(gen.fIsLaser)
+ ,fPrereadSCList(gen.fPrereadSCList)
{
//
gRandom->SetSeed();
//________________________________________________________________
AliToyMCEventGenerator::~AliToyMCEventGenerator()
{
- delete fTPCCorrection;
+ if (HasSCList() &&!fPrereadSCList) delete fTPCCorrection;
+ delete fSCList;
}
//________________________________________________________________
// this should be called after the tree was connected
//
+ if (HasSCList()) {
+ InitSpaceChargeList();
+ return;
+ }
+
AliInfo(Form("Using space charge map file: '%s'",fCorrectionFile.Data()));
TString corrName("map");
fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
}
}
+
+//________________________________________________________________
+void AliToyMCEventGenerator::InitSpaceChargeList()
+{
+ //
+ // init space charge conditions from a list of files
+ // this should be called after the tree was connected
+ //
+
+ std::ifstream file(fSCListFile.Data());
+ TString list;
+ list.ReadFile(file);
+
+ TObjArray *arr=list.Tokenize("\n");
+ if (!arr->GetEntriesFast()) {
+ delete arr;
+ AliFatal(Form("No SC file initialised. SC list '%s' seems empty",fSCListFile.Data()));
+ return;
+ }
+
+ // it is assumed that in case of an input list
+ // fCorrectionFile contains the name of the average correction
+ // to be then used in the reconstruction
+ if (fOutTree){
+ if (fCorrectionFile.IsNull()) {
+ AliFatal("List of SC files set, but no average map is specified. Use 'SetSpaceChargeFile' to do so");
+ return;
+ }
+ AliInfo("Attaching average space charge map file name to the tree");
+ fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
+ }
+
+
+ // case of non preread
+ // store the names of the files
+ if (!fPrereadSCList) {
+ if (fSCList) delete fSCList;
+ fSCList=arr;
+ return;
+ }
+
+ // case of preread
+ // load all SC files and set them to the list
+ if (!fSCList) fSCList=new TObjArray;
+ fSCList->SetOwner();
+ fSCList->Delete();
+
+ for (Int_t ifile=0; ifile<arr->GetEntriesFast(); ++ifile) {
+ TFile f(arr->At(ifile)->GetName());
+ if (!f.IsOpen()||f.IsZombie()) continue;
+ AliTPCSpaceCharge3D *sc=(AliTPCSpaceCharge3D*)f.Get("map");
+ if (!sc) continue;
+ sc->SetName(f.GetName());
+ fSCList->Add(sc);
+ }
+}
+
+//________________________________________________________________
+void AliToyMCEventGenerator::IterateSC(Int_t ipos)
+{
+ //
+ // In case a list of SC files is set iterate over them
+ //
+
+ if (!HasSCList()) return;
+
+ if (ipos<0) {
+ if (fOutTree) ipos=fOutTree->GetEntries();
+ else ipos=0;
+ }
+
+ TObject *sc=fSCList->At(ipos%fSCList->GetEntriesFast());
+ AliInfo(Form("Event: %d - SC: %s",ipos,sc->GetName()));
+ // case SC files have been preread
+ if (fPrereadSCList) {
+ fTPCCorrection=(AliTPCCorrection*)sc;
+ return;
+ }
+
+ // case no preread was done
+ TString &file=((TObjString*)sc)->String();
+ delete fTPCCorrection;
+ fTPCCorrection=0x0;
+
+ TFile f(file);
+ if (!f.IsOpen()||f.IsZombie()) {
+ AliFatal(Form("Could not open SC file '%s'",file.Data()));
+ return;
+ }
+
+ fTPCCorrection=(AliTPCSpaceCharge3D*)f.Get("map");
+ if (!fTPCCorrection) {
+ AliFatal(Form("Could not read SC map from SC file '%s'",file.Data()));
+ return;
+ }
+
+}
class TFile;
class TTree;
+class TObjArray;
class AliTPCParam;
class AliTPCCorrection;
void SetIsLaser(Bool_t use) { fIsLaser=use; }
Bool_t GetIsLaser() const { return fIsLaser; }
+
+ void SetSCListFile(const char* file) { fSCListFile=file; }
+ const char* GetSCListFile() const { return fSCListFile.Data(); }
+ void SetPrereadSCList(Bool_t b) { fPrereadSCList=b; }
+ Bool_t GetPrereadSCList() const { return fPrereadSCList; }
+ Bool_t HasSCList() const { return fSCList!=0x0; }
protected:
- AliTPCParam *fTPCParam;
- AliToyMCEvent *fEvent;
+ AliTPCParam *fTPCParam; //! TPC params
+ AliToyMCEvent *fEvent; //! Toy event
Bool_t ConnectOutputFile();
Bool_t CloseOutputFile();
void FillTree();
+ void IterateSC(Int_t ipos=-1);
UInt_t fCurrentTrack; // unique track id within the current event generation
private:
AliToyMCEventGenerator& operator= (const AliToyMCEventGenerator& );
- AliTPCCorrection *fTPCCorrection;
+ AliTPCCorrection *fTPCCorrection; //! distortion correction
+
+ TObjArray *fSCList; //! list with
+ TString fSCListFile; // file with a list of space charge files
+ TString fCorrectionFile; // name of a sinfle SC file
+ TString fOutputFileName; // name of the output file
+ TFile *fOutFile; //! output file
+ TTree *fOutTree; //! output tree
- TString fCorrectionFile;
- TString fOutputFileName;
- TFile *fOutFile;
- TTree *fOutTree;
+ Bool_t fUseStepCorrection; // use integralDz method?
+ Bool_t fUseMaterialBudget; // use material budget in tracking?
+ Bool_t fIsLaser; // is a laser event?
+ Bool_t fPrereadSCList; // preread all SC files from the SC list
- Bool_t fUseStepCorrection;
- Bool_t fUseMaterialBudget;
- Bool_t fIsLaser;
+ void InitSpaceChargeList();
ClassDef(AliToyMCEventGenerator, 1)
-
+
};
#endif
testExtrapolation();
*/
+#include "TStyle.h"
+#include "TCut.h"
+#include "TCanvas.h"
#include "TMath.h"
#include "TVectorD.h"
#include "TLinearFitter.h"
delete pcstream;
}
delete f;
+
+ gStyle->SetOptTitle(0);
f = TFile::Open("testExtrapolationErr.root","update");
TTree * tree = (TTree*)f->Get("extrapol");
//
TGraphErrors * grITSTRD0 = TStatToolkit::MakeGraphErrors(tree,"10*vecITSTRDErr0.fElements:vecR.fElements","",21,2,0);
//
grITS0->GetXaxis()->SetTitle("radius (cm)");
- grITS0->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
+ grITS0->GetYaxis()->SetTitle("#sigma_{r#varphi} (mm)");
grITS0->Draw("ap");
grITSTRD0->Draw("p");
- TLegend * legend = new TLegend(0.11,0.6,0.5,0.89,"Track residuals at TPC (q/p_{t}=0)");
+ TLegend * legend = new TLegend(0.11,0.65,0.55,0.89,"Track residuals at TPC (q/p_{t}=0)");
legend->AddEntry(grITS0,"ITS extrapolation","p");
legend->AddEntry(grITSTRD0,"ITS-TRD interpolation","p");
+ legend->SetFillColor(10);
legend->Draw();
//
//