[u/mrichter/AliRoot.git] / ITS / UPGRADE / testDetKMC.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TCanvas.h>
#include <TString.h>
#include <TROOT.h>
#include <TSystem.h>
#include <TFile.h>
#include <TObjArray.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include "KMCDetector.h"
#endif

TObjArray* CreateTrackConditions(const char* nm);
TObjArray  summary;


TObject* gDet=0;

// ptBins to generate
//double pts[] = {0.15,0.175,0.2,0.225,0.25,0.3,0.4,0.5,0.8,1.1,1.6,2,4,8,13,20};
double pts[] = {0.2};
//double pts[] = {0.2, 1., 10.};

void testDetKMC(int nev=1000,  // n events to generate per bin
		int setVer=0,   // optional version to pass for detector building
		double mass=0.14,  // particle mass to generate
		double eta=0.45,   // eta to test
		Bool_t aliceNew=kFALSE,
		Bool_t tpc=kTRUE, Double_t eff=0.95,  // detector settings
		Double_t vtxConstr=-1., // use vertex constraint in the fit
		const char* sumOut="trsumDef.root" // output file name
		)
{
  if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
  //
  // Detector, tracking initialization >>>
  //
  if (vtxConstr>0) KMCDetector::SetVtxConstraint(vtxConstr,vtxConstr);
  KMCLayer::SetDefEff(eff);
  //
  KMCDetector *itsdet = new KMCDetector((char*)"ALICE",(char*)"ITS");
  gDet = itsdet;
  KMCDetector &its = *itsdet;
  its.SetAvgRapidity(eta);
  if (aliceNew)     its.MakeAliceAllNew(tpc,1, setVer); // its sa
  else              its.MakeAliceCurrent(tpc,0); // its sa  
  //its.MakeAliceCurrent(0,1); // with tpc
  its.SetMaxSeedToPropagate(300);
  its.SetUseBackground();
  //
  // min hits per track to validate it
  its.RequireMinITSHits( its.GetNActiveITSLayers() - 3 );//(int)TMath::Max(4.0,(3./4.)*its.GetNActiveITSLayers()));
  printf("Min number of hits requested: %d\n",its.GetMinITSHits());
  //
  // max cluster-track chi2 
  its.RequireMaxChi2Cl(16.);
  // max chi2/NDF of the Kalman fit
  its.RequireMaxNormChi2NDF(5.);
  // penalty to chi2 from missing hit
  KMCProbe::SetMissingHitPenalty(6);
  //
  // Detector, tracking initialization <<<

  its.InitMCWatchHistos();
  summary.Clear();
  //
  int npt = sizeof(pts)/sizeof(double);
  printf("NPt bins = %d\n",npt);
  double ptMin = pts[0];
  double ptMax = pts[npt-1];
  //
  for (int ip=npt;ip--;) {
    //
    int nevR = nev;
    double pt = pts[ip];//1./(ptminI+dpt*ip);  // ptmin+ip*(ptmax-ptmin)/npt;

    double chi2Cl = 16.;
    if (npt>1) chi2Cl -= (16.-9.)*(1./pt-1./ptMax)/(1./ptMin-1./ptMax);
    //    if (pt<0.3) nevR = int(0.25*nevR);
    //    else if (pt<0.5) nevR = int(0.7*nevR);
    printf("Doing pt(%d)=%.3f for mass %.3f, %d ev | chi2Cl=%.2f\n",ip,pt,mass,nevR,chi2Cl);
    TString smn = Form("pt%d",ip);
    TObjArray* trSum = CreateTrackConditions(smn.Data());
    if (trSum) summary.AddLast(trSum);
    //
    its.SolveSingleTrack(mass, pt, eta, trSum, nevR);
  }
  //
  TString smout = sumOut;
  if (!smout.IsNull()) {
    gSystem->ExpandPathName(smout);
    TFile* fl = TFile::Open(smout.Data(),"recreate");
    if (!fl) {printf("Failed to open %s file\n",smout.Data()); return;}
    fl->WriteObject(&summary,"trSum","kSingleKey");
    fl->Close();
    delete fl;
  }
}

//_____________________________________________________________
TObjArray* CreateTrackConditions(const char* nm)
{
  // create set of track criteria to extract the summary
  KMCTrackSummary *sm = 0;
  TObjArray* arr = new TObjArray();
  //
  int nlr = KMCProbe::GetNITSLayers();
  //-------------------- put here user code -------------------
  // summary for ideal correct tracks
  sm = new KMCTrackSummary("corrAll","corrAll",nlr);
  sm->SetMinMaxClITS(nlr);
  sm->SetMinMaxClITSFake(0,0);
  sm->SetNamePrefix(nm);
  sm->AddPatternITSCorr( sm->Bits(1,0));    // the innermost 2 layers must have correct cluster
  sm->AddPatternITSCorr( sm->Bits(0,1));
  arr->AddLast(sm);
  //
  sm = new KMCTrackSummary("any","any",nlr);
  sm->SetMinMaxClITS(0);
  sm->SetMinMaxClITSFake(0,999);
  sm->SetNamePrefix(nm);
  arr->AddLast(sm);
  //
  // summary for good correct tracks
  sm = new KMCTrackSummary("corr3hSPD","corr3hSPD",nlr);
  sm->SetMinMaxClITS(3);
  sm->SetMinMaxClITSFake(0,0);
  sm->SetNamePrefix(nm);
  arr->AddLast(sm);
  //
  // summary for good correct tracks
  sm = new KMCTrackSummary("corr4h","corr4h",nlr);
  sm->SetMinMaxClITS(4);
  sm->SetMinMaxClITSFake(0,0);
  sm->SetNamePrefix(nm);
  arr->AddLast(sm);
  //
  // summary for good correct tracks with enough points at inner single layers
  sm = new KMCTrackSummary("corr4hs2in","corr4hs2in",nlr);
  sm->SetMinMaxClITS(4);
  sm->SetMinMaxClITSFake(0,0);
  sm->SetNamePrefix(nm);
  sm->AddPatternITSCorr( sm->Bits(1,0));    // the innermost 2 layers must have correct cluster
  sm->AddPatternITSCorr( sm->Bits(0,1));   
  arr->AddLast(sm);
  //
  // summary for good correct tracks with enough points at inner double layers
  sm = new KMCTrackSummary("corr4hd2in","corr5hd2in",nlr);
  sm->SetMinMaxClITS(4);
  sm->SetMinMaxClITSFake(0,0);
  sm->SetNamePrefix(nm);
  sm->AddPatternITSCorr( sm->Bits(1));    // the innermost 2 double layers must have correct cluster
  sm->AddPatternITSCorr( sm->Bits(0,1,1));   
  arr->AddLast(sm);
  //
  //
  // summary for golden tracks with AllNew RS setup: good pattern for offset and pt measurement
  sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
  sm->SetMinMaxClITS(4);
  sm->SetMinMaxClITSFake(0,0);
  sm->SetNamePrefix(nm);
  sm->AddPatternITSCorr( sm->Bits(1,1));    // the innermost 2 double layers must have correct cluster
  sm->AddPatternITSCorr( sm->Bits(0,0,1)); 
  //
  sm->AddPatternITSCorr( sm->Bits(0,0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
  //
  sm->AddPatternITSCorr( sm->Bits(0,0,0, 0,0, 1,1)); // measurement at large radius
  //
  arr->AddLast(sm);
  //
  // summary for golden tracks with curr setup: good pattern for offset and pt measurement
  sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
  sm->SetMinMaxClITS(4);
  sm->SetMinMaxClITSFake(0,0);
  sm->SetNamePrefix(nm);
  sm->AddPatternITSCorr( sm->Bits(1,0));    // the innermost 2 double layers must have correct cluster
  sm->AddPatternITSCorr( sm->Bits(0,1)); 
  //
  sm->AddPatternITSCorr( sm->Bits(0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
  //
  sm->AddPatternITSCorr( sm->Bits(0,0, 0,0, 1,1)); // measurement at large radius
  //
  arr->AddLast(sm);
  //
  // summary for tracks with at least 1 fake cluster
  sm = new KMCTrackSummary("fakeAny","fakeAny",nlr);
  sm->SetMinMaxClITS(4);
  sm->SetMinMaxClITSFake(1);
  sm->SetNamePrefix(nm);
  arr->AddLast(sm);
  //  
  // // summary for tracks with at least 1 fake cluster (but still acceptable?)
  // sm = new KMCTrackSummary("corrF","corrF",nlr);
  // sm->SetMinMaxClITS(4);
  // sm->SetMinMaxClITSCorr(3);              // at least 3 correct clusters
  // sm->SetMinMaxClITSFake(1);              // at least 1 fake clusters
  // //  /*
  // sm->AddPatternITSCorr( sm->Bits(1));    // the innermost 2 layers must have correct cluster
  // sm->AddPatternITSCorr( sm->Bits(0,1,1));   
  // //  */
  // //  sm->AddPatternITSCorr( sm->Bits(1,1));    // 
  // // sm->AddPatternITSCorr( sm->Bits(0,0,1,1,1,1));   

  // //  sm->AddPatternITSFakeExcl( sm->Bits(1,1,1)); // no fakes in innermost 3 layers
  // sm->SetNamePrefix(nm);
  // arr->AddLast(sm);
  //  
  //-----------------------------------------------------------
  return arr->GetEntries()>0 ? arr : 0;
}
Commit	Line	Data
c71b333a	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <TCanvas.h>
	3	#include <TString.h>
	4	#include <TROOT.h>
	5	#include <TSystem.h>
	6	#include <TFile.h>
	7	#include <TObjArray.h>
	8	#include <TH1.h>
	9	#include <TH2.h>
	10	#include <TF1.h>
	11	#include "KMCDetector.h"
	12	#endif
	13
	14	TObjArray* CreateTrackConditions(const char* nm);
	15	TObjArray summary;
	16
	17
	18	TObject* gDet=0;
	19
	20	// ptBins to generate
4f976c08	21	//double pts[] = {0.15,0.175,0.2,0.225,0.25,0.3,0.4,0.5,0.8,1.1,1.6,2,4,8,13,20};
4f976c08	22	double pts[] = {0.2};
c71b333a	23	//double pts[] = {0.2, 1., 10.};
c71b333a	24
4f976c08	25	void testDetKMC(int nev=1000, // n events to generate per bin
c71b333a	26	int setVer=0, // optional version to pass for detector building
	27	double mass=0.14, // particle mass to generate
	28	double eta=0.45, // eta to test
4f976c08	29	Bool_t aliceNew=kFALSE,
4f976c08	30	Bool_t tpc=kTRUE, Double_t eff=0.95, // detector settings
c71b333a	31	Double_t vtxConstr=-1., // use vertex constraint in the fit
	32	const char* sumOut="trsumDef.root" // output file name
	33	)
	34	{
	35	if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
	36	//
	37	// Detector, tracking initialization >>>
	38	//
	39	if (vtxConstr>0) KMCDetector::SetVtxConstraint(vtxConstr,vtxConstr);
	40	KMCLayer::SetDefEff(eff);
	41	//
	42	KMCDetector itsdet = new KMCDetector((char)"ALICE",(char*)"ITS");
	43	gDet = itsdet;
	44	KMCDetector &its = *itsdet;
	45	its.SetAvgRapidity(eta);
	46	if (aliceNew) its.MakeAliceAllNew(tpc,1, setVer); // its sa
	47	else its.MakeAliceCurrent(tpc,0); // its sa
	48	//its.MakeAliceCurrent(0,1); // with tpc
	49	its.SetMaxSeedToPropagate(300);
	50	its.SetUseBackground();
	51	//
	52	// min hits per track to validate it
4f976c08	53	its.RequireMinITSHits( its.GetNActiveITSLayers() - 3 );//(int)TMath::Max(4.0,(3./4.)*its.GetNActiveITSLayers()));
c71b333a	54	printf("Min number of hits requested: %d\n",its.GetMinITSHits());
	55	//
	56	// max cluster-track chi2
	57	its.RequireMaxChi2Cl(16.);
	58	// max chi2/NDF of the Kalman fit
	59	its.RequireMaxNormChi2NDF(5.);
	60	// penalty to chi2 from missing hit
	61	KMCProbe::SetMissingHitPenalty(6);
	62	//
	63	// Detector, tracking initialization <<<
	64
	65	its.InitMCWatchHistos();
	66	summary.Clear();
	67	//
	68	int npt = sizeof(pts)/sizeof(double);
4f976c08	69	printf("NPt bins = %d\n",npt);
c71b333a	70	double ptMin = pts[0];
	71	double ptMax = pts[npt-1];
	72	//
	73	for (int ip=npt;ip--;) {
	74	//
	75	int nevR = nev;
	76	double pt = pts[ip];//1./(ptminI+dptip); // ptmin+ip(ptmax-ptmin)/npt;
	77
4f976c08	78	double chi2Cl = 16.;
4f976c08	79	if (npt>1) chi2Cl -= (16.-9.)*(1./pt-1./ptMax)/(1./ptMin-1./ptMax);
c71b333a	80	// if (pt<0.3) nevR = int(0.25*nevR);
	81	// else if (pt<0.5) nevR = int(0.7*nevR);
	82	printf("Doing pt(%d)=%.3f for mass %.3f, %d ev \| chi2Cl=%.2f\n",ip,pt,mass,nevR,chi2Cl);
	83	TString smn = Form("pt%d",ip);
	84	TObjArray* trSum = CreateTrackConditions(smn.Data());
	85	if (trSum) summary.AddLast(trSum);
	86	//
	87	its.SolveSingleTrack(mass, pt, eta, trSum, nevR);
	88	}
	89	//
	90	TString smout = sumOut;
	91	if (!smout.IsNull()) {
	92	gSystem->ExpandPathName(smout);
	93	TFile* fl = TFile::Open(smout.Data(),"recreate");
	94	if (!fl) {printf("Failed to open %s file\n",smout.Data()); return;}
	95	fl->WriteObject(&summary,"trSum","kSingleKey");
	96	fl->Close();
	97	delete fl;
	98	}
	99	}
	100
	101	//_____________________________________________________________
	102	TObjArray* CreateTrackConditions(const char* nm)
	103	{
	104	// create set of track criteria to extract the summary
	105	KMCTrackSummary *sm = 0;
	106	TObjArray* arr = new TObjArray();
	107	//
	108	int nlr = KMCProbe::GetNITSLayers();
	109	//-------------------- put here user code -------------------
	110	// summary for ideal correct tracks
	111	sm = new KMCTrackSummary("corrAll","corrAll",nlr);
	112	sm->SetMinMaxClITS(nlr);
	113	sm->SetMinMaxClITSFake(0,0);
	114	sm->SetNamePrefix(nm);
4f976c08	115	sm->AddPatternITSCorr( sm->Bits(1,0)); // the innermost 2 layers must have correct cluster
	116	sm->AddPatternITSCorr( sm->Bits(0,1));
	117	arr->AddLast(sm);
	118	//
	119	sm = new KMCTrackSummary("any","any",nlr);
	120	sm->SetMinMaxClITS(0);
	121	sm->SetMinMaxClITSFake(0,999);
	122	sm->SetNamePrefix(nm);
	123	arr->AddLast(sm);
	124	//
	125	// summary for good correct tracks
	126	sm = new KMCTrackSummary("corr3hSPD","corr3hSPD",nlr);
	127	sm->SetMinMaxClITS(3);
	128	sm->SetMinMaxClITSFake(0,0);
	129	sm->SetNamePrefix(nm);
c71b333a	130	arr->AddLast(sm);
	131	//
	132	// summary for good correct tracks
	133	sm = new KMCTrackSummary("corr4h","corr4h",nlr);
	134	sm->SetMinMaxClITS(4);
	135	sm->SetMinMaxClITSFake(0,0);
	136	sm->SetNamePrefix(nm);
	137	arr->AddLast(sm);
	138	//
	139	// summary for good correct tracks with enough points at inner single layers
	140	sm = new KMCTrackSummary("corr4hs2in","corr4hs2in",nlr);
	141	sm->SetMinMaxClITS(4);
	142	sm->SetMinMaxClITSFake(0,0);
	143	sm->SetNamePrefix(nm);
	144	sm->AddPatternITSCorr( sm->Bits(1,0)); // the innermost 2 layers must have correct cluster
	145	sm->AddPatternITSCorr( sm->Bits(0,1));
	146	arr->AddLast(sm);
	147	//
	148	// summary for good correct tracks with enough points at inner double layers
	149	sm = new KMCTrackSummary("corr4hd2in","corr5hd2in",nlr);
	150	sm->SetMinMaxClITS(4);
	151	sm->SetMinMaxClITSFake(0,0);
	152	sm->SetNamePrefix(nm);
	153	sm->AddPatternITSCorr( sm->Bits(1)); // the innermost 2 double layers must have correct cluster
	154	sm->AddPatternITSCorr( sm->Bits(0,1,1));
	155	arr->AddLast(sm);
	156	//
	157	//
	158	// summary for golden tracks with AllNew RS setup: good pattern for offset and pt measurement
	159	sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
	160	sm->SetMinMaxClITS(4);
	161	sm->SetMinMaxClITSFake(0,0);
	162	sm->SetNamePrefix(nm);
	163	sm->AddPatternITSCorr( sm->Bits(1,1)); // the innermost 2 double layers must have correct cluster
	164	sm->AddPatternITSCorr( sm->Bits(0,0,1));
	165	//
	166	sm->AddPatternITSCorr( sm->Bits(0,0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
	167	//
	168	sm->AddPatternITSCorr( sm->Bits(0,0,0, 0,0, 1,1)); // measurement at large radius
	169	//
	170	arr->AddLast(sm);
	171	//
	172	// summary for golden tracks with curr setup: good pattern for offset and pt measurement
	173	sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
	174	sm->SetMinMaxClITS(4);
	175	sm->SetMinMaxClITSFake(0,0);
	176	sm->SetNamePrefix(nm);
	177	sm->AddPatternITSCorr( sm->Bits(1,0)); // the innermost 2 double layers must have correct cluster
	178	sm->AddPatternITSCorr( sm->Bits(0,1));
	179	//
	180	sm->AddPatternITSCorr( sm->Bits(0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
	181	//
	182	sm->AddPatternITSCorr( sm->Bits(0,0, 0,0, 1,1)); // measurement at large radius
	183	//
	184	arr->AddLast(sm);
	185	//
	186	// summary for tracks with at least 1 fake cluster
	187	sm = new KMCTrackSummary("fakeAny","fakeAny",nlr);
	188	sm->SetMinMaxClITS(4);
	189	sm->SetMinMaxClITSFake(1);
	190	sm->SetNamePrefix(nm);
	191	arr->AddLast(sm);
	192	//
	193	// // summary for tracks with at least 1 fake cluster (but still acceptable?)
194	// sm = new KMCTrackSummary("corrF","corrF",nlr);
195	// sm->SetMinMaxClITS(4);
196	// sm->SetMinMaxClITSCorr(3); // at least 3 correct clusters
197	// sm->SetMinMaxClITSFake(1); // at least 1 fake clusters
198	// // /*
199	// sm->AddPatternITSCorr( sm->Bits(1)); // the innermost 2 layers must have correct cluster
200	// sm->AddPatternITSCorr( sm->Bits(0,1,1));
201	// // */
202	// // sm->AddPatternITSCorr( sm->Bits(1,1)); //
203	// // sm->AddPatternITSCorr( sm->Bits(0,0,1,1,1,1));
204
205	// // sm->AddPatternITSFakeExcl( sm->Bits(1,1,1)); // no fakes in innermost 3 layers
206	// sm->SetNamePrefix(nm);
207	// arr->AddLast(sm);
208	//
209	//-----------------------------------------------------------
210	return arr->GetEntries()>0 ? arr : 0;
211	}