[u/mrichter/AliRoot.git] / MUON / mapping / macros / testSegmentation.C

// $Id$
// $MpId: testSegmentation.C,v 1.1 2005/09/19 19:02:53 ivana Exp $

#include "AliRun.h"
#include "AliMUON.h"
#include "AliMUONChamber.h"
#include "AliMUONGeometrySegmentation.h"
#include "AliMUONVGeometryDESegmentation.h"
#include "AliMpPlaneType.h"
#include "AliMUONSt345SlatSegmentationV2.h"
//#include "AliMpSt345Reader.h"
#include "AliMpSegmentationManager.h"
#include "AliMpSlat.h"
#include "AliMpPCB.h"

#include <iostream>
#include <map>
#include <set>
#include <string>
#include <sstream>

std::map<int,std::string> detElemIdToSlatTypeMap;

//_____________________________________________________________________________
bool equal(double a, double b, double epsilon=1e-5)
{
  double d = fabs(a-b)/b;
  if ( d < epsilon ) return true;
  return false;
}

//_____________________________________________________________________________
bool error(const char* msg)
{
  std::cout << msg << std::endl;
  return false;
}

//_____________________________________________________________________________
bool readDetElemIdToSlatType(const char* file = "/Users/aphecetc/Workspace/aliroot/MUON/mapping/data/station345/DetElemIdToSlatType.dat")
{
  std::ifstream in(file);
  if (!in.good()) return false;

  char line[80];

  while ( in.getline(line,80) )
    {
      if ( !isdigit(line[0]) ) continue;
      std::string sline(line);
      
      std::string::size_type pos = sline.find_first_of(' ');
      int detelemid = std::atoi(sline.substr(0,pos).c_str());
      detElemIdToSlatTypeMap[detelemid] = sline.substr(pos+1);
    }

  in.close();

  return true;
}

//_____________________________________________________________________________
const char* SlatType(int detelemid)
{
  if ( detElemIdToSlatTypeMap.empty() ) readDetElemIdToSlatType();

  std::map<int,std::string>::const_iterator it = 
    detElemIdToSlatTypeMap.find(detelemid);

  if ( it != detElemIdToSlatTypeMap.end() )
    {
      return it->second.c_str();
    }
  else
    {
      return 0;
    }
}

//_____________________________________________________________________________
bool testAllDetelem()
{
  // test we can read the map file (detelemid <-> slat type)
  bool ok = readDetElemIdToSlatType();
  if (!ok)
    {
      return false;
    }

  // now loops over all detelemid and check that:
  // 1. We can actually read the slat files, both bending and non-bending
  // 2. Bending and non-bending slats do have the same x-size

  std::map<int,std::string>::const_iterator it;
  std::set<std::string> slattypes;
  std::map<std::string,std::string> problems;

  for ( it = detElemIdToSlatTypeMap.begin(); it != detElemIdToSlatTypeMap.end();
	++it )
    {
      slattypes.insert(it->second);
      std::cout << it->first << " : " << it->second <<  " ";
      AliMpSlat* b = static_cast<AliMpSlat*>
	(AliMpSegmentationManager::Segmentation(it->first,AliMp::kBendingPlane));
      AliMpSlat* nb = static_cast<AliMpSlat*>
	(AliMpSegmentationManager::Segmentation(it->first,AliMp::kNonBendingPlane));
      std::cout << " B : " << b << " NB : " << nb;
      Double_t bx = 0.0;
      Double_t nbx = 0.0;
      std::ostringstream pb;

      if (b) 
	{
	  bx = b->DX()*2;
	}
      else
	{
	  std::cout << " Missing BENDING. ";
	}
     if (nb) 
	{
	  nbx = nb->DX()*2;
	}
     else
       {
	 std::cout << " Missing NON-BENDING.";
       }
     std::cout << std::endl;
     if (!b || !nb) 
       {
	 return false;
       }
     if ( bx != nbx )
       {
	 // Find which pcb(s) have a different size
	 if ( b->GetSize() != nb->GetSize() )
	   {
	     std::cout << "Not the same number of PCBs !" << std::endl;
	     return false;
	   }

	 pb << " DIFFERENT SIZES ! Bending= " << bx
	    << " Non-Bending= " << nbx
	    << " delta = " << bx - nbx;

	 for ( size_t i = 0; i < b->GetSize(); ++ i )
	   {
	     Double_t delta = b->GetPCB(i)->DX() - nb->GetPCB(i)->DX();
	     if ( fabs(delta*2) > 1E-3 )
	       {
		 pb << " DELTA(" << i << ")=" << delta*2;
	       }
	   }	 
	 pb << std::ends;
	 problems[it->second] = pb.str();
       }
    }

  //
  std::cout << "Number of detelem (per plane): " 
	    << detElemIdToSlatTypeMap.size()
	    << std::endl;
  std::cout << "Number of slat types : " << slattypes.size()
	    << std::endl;
  std::cout << "(Potential) Problems:" << std::endl;

  std::map<std::string,std::string>::const_iterator sit;

  for ( sit = problems.begin(); sit != problems.end(); ++sit ) 
    {
      std::cout << sit->first << " : " << sit->second << std::endl;
    }

  return ok;
}

//_____________________________________________________________________________
bool getSegmentation(int detElemId, 
		     AliMpPlaneType bendingornot,
		     AliMUONVGeometryDESegmentation*& oldSegm,
		     AliMUONVGeometryDESegmentation*& newSegm)
{
  AliMUON* muon = static_cast<AliMUON*>(gAlice->GetModule("MUON"));
  if (!muon) 
    {
      gAlice->Init("${ALICE_ROOT}/MUON/Config.C");
      muon = static_cast<AliMUON*>(gAlice->GetModule("MUON"));
      if (!muon) return error("Cannot get MUON !");
    }

  int ichamber = detElemId/100 - 1;
  //  std::cout << "muon->Chamber(" << ichamber << ")" << std::endl;
  AliMUONChamber& chamber = muon->Chamber(ichamber);

  int icathode = 1;
  if ( bendingornot == AliMp::kNonBendingPlane ) icathode = 2;

  AliMUONGeometrySegmentation* gs = chamber.SegmentationModel2(icathode);
  if (!gs) return error(Form("Cannot get cathode %d",icathode));

  oldSegm = const_cast<AliMUONVGeometryDESegmentation*>(gs->GetDESegmentation(detElemId));
  if (!oldSegm) 
    {
      return error(Form("Cannot get segmentation for detElemId",detElemId));
    }

  newSegm =
    new AliMUONSt345SlatSegmentationV2(detElemId,bendingornot);

  return true;
}

//_____________________________________________________________________________
Int_t countPads(AliMUONVGeometryDESegmentation* seg)
{
  if (!seg) return 0;

  Int_t npads = 0;

  for ( Int_t ix = 1; ix <= seg->Npx(); ++ix )
    {
      for ( Int_t iy = 1; iy <= seg->Npy(); ++iy )
	{
	  if ( seg->HasPad(ix,iy) ) ++npads;
	}
    }
  return npads;
}

//_____________________________________________________________________________
bool testIC(Int_t d1, Int_t d2, AliMpPlaneType bendingornot)
{
  readDetElemIdToSlatType();
  AliMUONVGeometryDESegmentation* o;
  AliMUONVGeometryDESegmentation* n;

  Int_t ntested = 0;

  for ( Int_t d = d1; d <= d2; ++d )
    {
      if ( detElemIdToSlatTypeMap.find(d) == detElemIdToSlatTypeMap.end() ) continue;
      std::cout << d << std::endl;
      bool ok = getSegmentation(d,bendingornot,o,n);
      if (!ok) return false;

      Int_t nx = std::max(o->Npx(),n->Npx());
      Int_t ny = std::max(o->Npy(),n->Npy());
      
      for ( Int_t ix = 1; ix <= nx; ++ix )
	{
	  for ( Int_t iy = 1; iy <= ny; ++iy )
	    {
	      float xn,yn,zn;
	      float xo,yo,zo;
	      if ( !o->HasPad(ix,iy) || !n->HasPad(ix,iy) ) continue;
	      o->GetPadC(ix,iy,xo,yo,zo);
	      n->GetPadC(ix,iy,xn,yn,zn);
	      ++ntested;
	      if ( !equal(xn,xo) || !equal(yn,yo) )
		{
		  printf("%4d (%4d,%4d) -> OLD (%e,%e) NEW (%e,%e) DELTA (%e,%e)\n",d,ix,iy,xo,yo,xn,yn,xn-xo,yn-yo);
		}
	      // Now tries to go back to ix,iy from positions
	      Int_t nix,niy;
	      Int_t oix,oiy;
	      o->GetPadI(xo,yo,zo,oix,oiy);
	      n->GetPadI(xn,yn,zn,nix,niy);
	      std::string msg;
	      if ( ix != oix || iy != oiy )
		{
		  msg += "OLD";
		}
	      if ( ix != nix || iy != niy )
		{
		  msg += "NEW";
		}
	      if ( !msg.empty() )
		{
		  printf("Circular test failed for %s : (%4d,%4d) -> OLD (%e,%e) NEW (%e,%e) -> OLD (%4d,%4d) NEW (%4d,%4d)\n",msg.c_str(),ix,iy,xo,yo,xn,yn,oix,oiy,nix,niy);		  
		  return false;
		}
	    }
	}      
      delete n;
    }
  std::cout << "Number of tested pads  = " << ntested << std::endl;

  return true;
}

//_____________________________________________________________________________
void countPads(Int_t d1, Int_t d2)
{
  readDetElemIdToSlatType();

  Int_t onpads = 0;
  Int_t nnpads = 0;
  AliMUONVGeometryDESegmentation* o;
  AliMUONVGeometryDESegmentation* n;

  AliMpPlaneType pt[] = { AliMp::kNonBendingPlane, AliMp::kBendingPlane };

  for ( int ipt = 0; ipt < 2; ++ipt )
    {
      std::cout << ( (pt[ipt] == AliMp::kNonBendingPlane )?"NonBending"
		     : "Bending" )
		<< std::endl;
      Int_t p_nnpads = 0;
      Int_t p_onpads = 0;
      for ( Int_t d = d1; d <= d2; ++d )
	{
	  if ( detElemIdToSlatTypeMap.find(d) == detElemIdToSlatTypeMap.end() )
	    continue; // skip non-existing detElemId
	  bool ok = getSegmentation(d,pt[ipt],o,n);
	  if (!ok) return;
	  Int_t nn = countPads(n);
	  Int_t oo = countPads(o);
	  p_nnpads += nn;
	  p_onpads += oo;
	  printf("%4d OLD %5d NEW %5d Delta %6d OLDIXMAX %3d NEWIXMAX %3d\n",d,oo,nn,nn-oo,(o?o->Npx():0),n->Npx());
	  delete n;
	}

      std::cout << "OLD:" << p_onpads << std::endl
		<< "NEW:" << p_nnpads << std::endl;
      nnpads += p_nnpads;
      onpads += p_onpads;
    }

  std::cout << "-----" << std::endl
	    << "OLD:" << onpads << std::endl
	    << "NEW:" << nnpads << std::endl;
}
Commit	Line	Data
dd3e6c2b	1	// $Id$
	2	// $MpId: testSegmentation.C,v 1.1 2005/09/19 19:02:53 ivana Exp $
	3
	4	#include "AliRun.h"
	5	#include "AliMUON.h"
	6	#include "AliMUONChamber.h"
	7	#include "AliMUONGeometrySegmentation.h"
	8	#include "AliMUONVGeometryDESegmentation.h"
	9	#include "AliMpPlaneType.h"
	10	#include "AliMUONSt345SlatSegmentationV2.h"
	11	//#include "AliMpSt345Reader.h"
	12	#include "AliMpSegmentationManager.h"
	13	#include "AliMpSlat.h"
	14	#include "AliMpPCB.h"
	15
	16	#include <iostream>
	17	#include <map>
	18	#include <set>
	19	#include <string>
	20	#include <sstream>
	21
	22	std::map<int,std::string> detElemIdToSlatTypeMap;
	23
	24	//_____________________________________________________________________________
	25	bool equal(double a, double b, double epsilon=1e-5)
	26	{
	27	double d = fabs(a-b)/b;
	28	if ( d < epsilon ) return true;
	29	return false;
	30	}
	31
	32	//_____________________________________________________________________________
	33	bool error(const char* msg)
	34	{
	35	std::cout << msg << std::endl;
	36	return false;
	37	}
	38
	39	//_____________________________________________________________________________
	40	bool readDetElemIdToSlatType(const char* file = "/Users/aphecetc/Workspace/aliroot/MUON/mapping/data/station345/DetElemIdToSlatType.dat")
	41	{
	42	std::ifstream in(file);
	43	if (!in.good()) return false;
	44
	45	char line[80];
	46
	47	while ( in.getline(line,80) )
	48	{
	49	if ( !isdigit(line[0]) ) continue;
	50	std::string sline(line);
	51
	52	std::string::size_type pos = sline.find_first_of(' ');
	53	int detelemid = std::atoi(sline.substr(0,pos).c_str());
	54	detElemIdToSlatTypeMap[detelemid] = sline.substr(pos+1);
	55	}
	56
	57	in.close();
	58
	59	return true;
	60	}
	61
	62	//_____________________________________________________________________________
	63	const char* SlatType(int detelemid)
	64	{
65	if ( detElemIdToSlatTypeMap.empty() ) readDetElemIdToSlatType();
66
67	std::map<int,std::string>::const_iterator it =
68	detElemIdToSlatTypeMap.find(detelemid);
69
70	if ( it != detElemIdToSlatTypeMap.end() )
71	{
72	return it->second.c_str();
73	}
74	else
75	{
76	return 0;
77	}
78	}
79
80	//_____________________________________________________________________________
81	bool testAllDetelem()
82	{
83	// test we can read the map file (detelemid <-> slat type)
84	bool ok = readDetElemIdToSlatType();
85	if (!ok)
86	{
87	return false;
88	}
89
90	// now loops over all detelemid and check that:
91	// 1. We can actually read the slat files, both bending and non-bending
92	// 2. Bending and non-bending slats do have the same x-size
93
94	std::map<int,std::string>::const_iterator it;
95	std::set<std::string> slattypes;
96	std::map<std::string,std::string> problems;
97
98	for ( it = detElemIdToSlatTypeMap.begin(); it != detElemIdToSlatTypeMap.end();
99	++it )
100	{
101	slattypes.insert(it->second);
102	std::cout << it->first << " : " << it->second << " ";
103	AliMpSlat* b = static_cast<AliMpSlat*>
cddd101e	104	(AliMpSegmentationManager::Segmentation(it->first,AliMp::kBendingPlane));
dd3e6c2b	105	AliMpSlat* nb = static_cast<AliMpSlat*>
cddd101e	106	(AliMpSegmentationManager::Segmentation(it->first,AliMp::kNonBendingPlane));
dd3e6c2b	107	std::cout << " B : " << b << " NB : " << nb;
	108	Double_t bx = 0.0;
	109	Double_t nbx = 0.0;
	110	std::ostringstream pb;
	111
	112	if (b)
	113	{
	114	bx = b->DX()*2;
	115	}
	116	else
	117	{
	118	std::cout << " Missing BENDING. ";
	119	}
	120	if (nb)
	121	{
	122	nbx = nb->DX()*2;
	123	}
	124	else
	125	{
	126	std::cout << " Missing NON-BENDING.";
	127	}
	128	std::cout << std::endl;
	129	if (!b \|\| !nb)
	130	{
	131	return false;
	132	}
	133	if ( bx != nbx )
	134	{
	135	// Find which pcb(s) have a different size
	136	if ( b->GetSize() != nb->GetSize() )
	137	{
	138	std::cout << "Not the same number of PCBs !" << std::endl;
	139	return false;
	140	}
	141
	142	pb << " DIFFERENT SIZES ! Bending= " << bx
	143	<< " Non-Bending= " << nbx
	144	<< " delta = " << bx - nbx;
	145
	146	for ( size_t i = 0; i < b->GetSize(); ++ i )
	147	{
	148	Double_t delta = b->GetPCB(i)->DX() - nb->GetPCB(i)->DX();
	149	if ( fabs(delta*2) > 1E-3 )
	150	{
	151	pb << " DELTA(" << i << ")=" << delta*2;
	152	}
	153	}
	154	pb << std::ends;
	155	problems[it->second] = pb.str();
	156	}
	157	}
	158
	159	//
	160	std::cout << "Number of detelem (per plane): "
	161	<< detElemIdToSlatTypeMap.size()
	162	<< std::endl;
	163	std::cout << "Number of slat types : " << slattypes.size()
	164	<< std::endl;
	165	std::cout << "(Potential) Problems:" << std::endl;
	166
	167	std::map<std::string,std::string>::const_iterator sit;
	168
	169	for ( sit = problems.begin(); sit != problems.end(); ++sit )
	170	{
171	std::cout << sit->first << " : " << sit->second << std::endl;
172	}
173
174	return ok;
175	}
176
177	//_____________________________________________________________________________
178	bool getSegmentation(int detElemId,
179	AliMpPlaneType bendingornot,
180	AliMUONVGeometryDESegmentation*& oldSegm,
181	AliMUONVGeometryDESegmentation*& newSegm)
182	{
183	AliMUON* muon = static_cast<AliMUON*>(gAlice->GetModule("MUON"));
184	if (!muon)
185	{
186	gAlice->Init("${ALICE_ROOT}/MUON/Config.C");
187	muon = static_cast<AliMUON*>(gAlice->GetModule("MUON"));
188	if (!muon) return error("Cannot get MUON !");
189	}
190
191	int ichamber = detElemId/100 - 1;
192	// std::cout << "muon->Chamber(" << ichamber << ")" << std::endl;
193	AliMUONChamber& chamber = muon->Chamber(ichamber);
194
195	int icathode = 1;
cddd101e	196	if ( bendingornot == AliMp::kNonBendingPlane ) icathode = 2;
dd3e6c2b	197
	198	AliMUONGeometrySegmentation* gs = chamber.SegmentationModel2(icathode);
	199	if (!gs) return error(Form("Cannot get cathode %d",icathode));
	200
	201	oldSegm = const_cast<AliMUONVGeometryDESegmentation*>(gs->GetDESegmentation(detElemId));
	202	if (!oldSegm)
	203	{
	204	return error(Form("Cannot get segmentation for detElemId",detElemId));
	205	}
	206
	207	newSegm =
	208	new AliMUONSt345SlatSegmentationV2(detElemId,bendingornot);
	209
	210	return true;
	211	}
	212
	213	//_____________________________________________________________________________
	214	Int_t countPads(AliMUONVGeometryDESegmentation* seg)
	215	{
	216	if (!seg) return 0;
	217
	218	Int_t npads = 0;
	219
	220	for ( Int_t ix = 1; ix <= seg->Npx(); ++ix )
	221	{
	222	for ( Int_t iy = 1; iy <= seg->Npy(); ++iy )
	223	{
	224	if ( seg->HasPad(ix,iy) ) ++npads;
	225	}
	226	}
	227	return npads;
	228	}
	229
	230	//_____________________________________________________________________________
	231	bool testIC(Int_t d1, Int_t d2, AliMpPlaneType bendingornot)
	232	{
	233	readDetElemIdToSlatType();
	234	AliMUONVGeometryDESegmentation* o;
	235	AliMUONVGeometryDESegmentation* n;
	236
	237	Int_t ntested = 0;
	238
	239	for ( Int_t d = d1; d <= d2; ++d )
	240	{
	241	if ( detElemIdToSlatTypeMap.find(d) == detElemIdToSlatTypeMap.end() ) continue;
	242	std::cout << d << std::endl;
	243	bool ok = getSegmentation(d,bendingornot,o,n);
	244	if (!ok) return false;
	245
	246	Int_t nx = std::max(o->Npx(),n->Npx());
	247	Int_t ny = std::max(o->Npy(),n->Npy());
	248
	249	for ( Int_t ix = 1; ix <= nx; ++ix )
	250	{
	251	for ( Int_t iy = 1; iy <= ny; ++iy )
	252	{
	253	float xn,yn,zn;
	254	float xo,yo,zo;
	255	if ( !o->HasPad(ix,iy) \|\| !n->HasPad(ix,iy) ) continue;
	256	o->GetPadC(ix,iy,xo,yo,zo);
	257	n->GetPadC(ix,iy,xn,yn,zn);
	258	++ntested;
	259	if ( !equal(xn,xo) \|\| !equal(yn,yo) )
	260	{
261	printf("%4d (%4d,%4d) -> OLD (%e,%e) NEW (%e,%e) DELTA (%e,%e)\n",d,ix,iy,xo,yo,xn,yn,xn-xo,yn-yo);
262	}
263	// Now tries to go back to ix,iy from positions
264	Int_t nix,niy;
265	Int_t oix,oiy;
266	o->GetPadI(xo,yo,zo,oix,oiy);
267	n->GetPadI(xn,yn,zn,nix,niy);
268	std::string msg;
269	if ( ix != oix \|\| iy != oiy )
270	{
271	msg += "OLD";
272	}
273	if ( ix != nix \|\| iy != niy )
274	{
275	msg += "NEW";
276	}
277	if ( !msg.empty() )
278	{
279	printf("Circular test failed for %s : (%4d,%4d) -> OLD (%e,%e) NEW (%e,%e) -> OLD (%4d,%4d) NEW (%4d,%4d)\n",msg.c_str(),ix,iy,xo,yo,xn,yn,oix,oiy,nix,niy);
280	return false;
281	}
282	}
283	}
284	delete n;
285	}
286	std::cout << "Number of tested pads = " << ntested << std::endl;
287
288	return true;
289	}
290
291	//_____________________________________________________________________________
292	void countPads(Int_t d1, Int_t d2)
293	{
294	readDetElemIdToSlatType();
295
296	Int_t onpads = 0;
297	Int_t nnpads = 0;
298	AliMUONVGeometryDESegmentation* o;
299	AliMUONVGeometryDESegmentation* n;
300
cddd101e	301	AliMpPlaneType pt[] = { AliMp::kNonBendingPlane, AliMp::kBendingPlane };
dd3e6c2b	302
	303	for ( int ipt = 0; ipt < 2; ++ipt )
	304	{
cddd101e	305	std::cout << ( (pt[ipt] == AliMp::kNonBendingPlane )?"NonBending"
dd3e6c2b	306	: "Bending" )
	307	<< std::endl;
	308	Int_t p_nnpads = 0;
	309	Int_t p_onpads = 0;
	310	for ( Int_t d = d1; d <= d2; ++d )
	311	{
	312	if ( detElemIdToSlatTypeMap.find(d) == detElemIdToSlatTypeMap.end() )
	313	continue; // skip non-existing detElemId
	314	bool ok = getSegmentation(d,pt[ipt],o,n);
	315	if (!ok) return;
	316	Int_t nn = countPads(n);
	317	Int_t oo = countPads(o);
	318	p_nnpads += nn;
	319	p_onpads += oo;
	320	printf("%4d OLD %5d NEW %5d Delta %6d OLDIXMAX %3d NEWIXMAX %3d\n",d,oo,nn,nn-oo,(o?o->Npx():0),n->Npx());
	321	delete n;
	322	}
	323
	324	std::cout << "OLD:" << p_onpads << std::endl
	325	<< "NEW:" << p_nnpads << std::endl;
	326	nnpads += p_nnpads;
	327	onpads += p_onpads;
	328	}
	329
	330	std::cout << "-----" << std::endl
	331	<< "OLD:" << onpads << std::endl
	332	<< "NEW:" << nnpads << std::endl;
	333	}