MUON/mapping/macros/testSegmentation.C

   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*>
 104         (AliMpSegmentationManager::Segmentation(it->first,AliMp::kBendingPlane));
 105       AliMpSlat* nb = static_cast<AliMpSlat*>
 106         (AliMpSegmentationManager::Segmentation(it->first,AliMp::kNonBendingPlane));
 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;
 196   if ( bendingornot == AliMp::kNonBendingPlane ) icathode = 2;
 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
 301   AliMpPlaneType pt[] = { AliMp::kNonBendingPlane, AliMp::kBendingPlane };
 302
 303   for ( int ipt = 0; ipt < 2; ++ipt )
 304     {
 305       std::cout << ( (pt[ipt] == AliMp::kNonBendingPlane )?"NonBending"
 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 }