MUON/mapping/macros/testSegmentation.C

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