MUON/AliMUONTriggerChamberEff.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /// \class AliMUONTriggerChamberEff
  19 /// implementation of the trigger chamber efficiency determination from
  20 /// data, and returns the
  21 /// efficiencyCells.dat with the calculated efficiencies
  22 ///
  23 /// \author Diego Stocco (Torino)
  24
  25 #include "AliMUONTriggerChamberEff.h"
  26 #include "AliMUONDigit.h"
  27 #include "AliMUONConstants.h"
  28 #include "AliMUONGlobalTrigger.h"
  29 #include "AliMUONGeometryTransformer.h"
  30 #include "AliMUONSegmentation.h"
  31 #include "AliMUON.h"
  32 #include "AliMUONData.h"
  33 #include "AliMUONTriggerTrack.h"
  34
  35 #include "AliMpVSegmentation.h"
  36 #include "AliMpSegmentation.h"
  37 #include "AliMpPad.h"
  38 #include "AliMpDEIterator.h"
  39 #include "AliMpPlaneType.h"
  40
  41 #include "AliRunLoader.h"
  42 #include "AliRun.h"
  43
  44 #include <Riostream.h>
  45 #include <TFile.h>
  46 #include <TH1F.h>
  47 #include <TMath.h>
  48
  49 /// \cond CLASSIMP
  50 ClassImp(AliMUONTriggerChamberEff)
  51 /// \endcond
  52
  53 //_____________________________________________________________________________
  54 AliMUONTriggerChamberEff::AliMUONTriggerChamberEff(const char* galiceFile,
  55                                                    Int_t firstEvent, Int_t lastEvent)
  56 : TObject(),
  57   fFirstEvent(firstEvent),
  58   fLastEvent(lastEvent),
  59   fFirstRun(-1),
  60   fLastRun(-1),
  61   fRunLoader(0x0),
  62   fData(0x0),
  63   fReproduceTrigResponse(kFALSE),
  64   fPrintInfo(kFALSE),
  65   fMUON(0x0),
  66   fDebugLevel(0),
  67   fGaliceDir(0x0)
  68 {
  69 /// Standard constructor
  70     SetGaliceFile(galiceFile);
  71     ResetArrays();
  72 }
  73
  74 //_____________________________________________________________________________
  75 AliMUONTriggerChamberEff::AliMUONTriggerChamberEff(Int_t firstRun, Int_t lastRun,
  76                                                    const char* galiceRunDir,
  77                                                    Int_t firstEvent, Int_t lastEvent)
  78 : TObject(),
  79   fFirstEvent(firstEvent),
  80   fLastEvent(lastEvent),
  81   fFirstRun(firstRun),
  82   fLastRun(lastRun),
  83   fRunLoader(0x0),
  84   fData(0x0),
  85   fReproduceTrigResponse(kFALSE),
  86   fPrintInfo(kFALSE),
  87   fMUON(0x0),
  88   fDebugLevel(0),
  89   fGaliceDir(galiceRunDir)
  90 {
  91 /// Standard constructor
  92     ResetArrays();
  93     delete gAlice;
  94 }
  95
  96 //_____________________________________________________________________________
  97 AliMUONTriggerChamberEff::~AliMUONTriggerChamberEff()
  98 {
  99 /// Destructor
 100     fRunLoader->UnloadAll();
 101     delete fRunLoader;
 102     delete fData;
 103 }
 104
 105 //_____________________________________________________________________________
 106 void AliMUONTriggerChamberEff::SetGaliceFile(const char *galiceFile)
 107 {
 108     //
 109     /// Opens the galice.root and loads tracks and digits.
 110     //
 111
 112     fRunLoader = AliRunLoader::Open(galiceFile,"MUONFolder","READ");
 113     if (!fRunLoader)
 114     {
 115         AliError(Form("Error opening %s file \n",galiceFile));
 116     }
 117     else
 118     {
 119         fRunLoader->LoadgAlice();
 120         gAlice = fRunLoader->GetAliRun();
 121         fMUON = (AliMUON*)gAlice->GetModule("MUON");
 122
 123         if(fLastEvent<=0 || fLastEvent>fRunLoader->GetNumberOfEvents())fLastEvent = fRunLoader->GetNumberOfEvents()-1;
 124         if(fFirstEvent<0)fFirstEvent=0;
 125
 126
 127         AliLoader* loader = fRunLoader->GetLoader("MUONLoader");
 128         if ( loader )
 129         {
 130             fData = new AliMUONData(loader,"MUON","MUON");
 131             loader->LoadTracks("READ");
 132             loader->LoadDigits("READ");
 133         }
 134         else
 135         {
 136             AliError(Form("Could get MUONLoader"));
 137         }
 138     }
 139 }
 140
 141 //_____________________________________________________________________________
 142 void AliMUONTriggerChamberEff::CleanGalice()
 143 {
 144     //
 145     /// Unload all loaded data
 146     //
 147
 148     fRunLoader->UnloadAll();
 149     delete fRunLoader;
 150     fRunLoader = 0;
 151 }
 152
 153 //_____________________________________________________________________________
 154 void AliMUONTriggerChamberEff::ResetArrays()
 155 {
 156     //
 157     /// Sets the data member counters to 0.
 158     //
 159
 160     for(Int_t ch=0; ch<fgkNchambers; ch++){
 161         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 162             fTrigger34[ch][cath] = 0;
 163             fTrigger44[cath] = 0;
 164             for(Int_t slat=0; slat<fgkNslats; slat++){
 165                 fInefficientSlat[ch][cath][slat] = 0;
 166                 fHitPerSlat[ch][cath][slat] = 0;
 167             }
 168         }
 169     }
 170 }
 171
 172
 173 //_____________________________________________________________________________
 174 void AliMUONTriggerChamberEff::InfoDigit()
 175 {
 176     //
 177     /// Prints information on digits (for debugging)
 178     //
 179
 180     AliMUONDigit * mDigit=0x0;
 181     Int_t firstTrigCh = AliMUONConstants::NTrackingCh();
 182     // Addressing
 183     Int_t nchambers = AliMUONConstants::NCh();
 184     fData->SetTreeAddress("D,GLT");
 185
 186     fData->GetDigits();
 187     // Loop on chambers
 188     for(Int_t ichamber=firstTrigCh; ichamber<nchambers; ichamber++) {
 189       TClonesArray* digits = fData->Digits(ichamber);
 190       digits->Sort();
 191       Int_t ndigits = (Int_t)digits->GetEntriesFast();
 192       for(Int_t idigit=0; idigit<ndigits; idigit++) {
 193           mDigit = (AliMUONDigit*)digits->At(idigit);
 194           mDigit->Print();
 195       } // end digit loop
 196     } // end chamber loop
 197     fData->ResetDigits();
 198 }
 199
 200
 201 //_____________________________________________________________________________
 202 Bool_t AliMUONTriggerChamberEff::PadMatchTrack(Float_t xPad, Float_t yPad, Float_t dpx, Float_t dpy,
 203                                                Float_t xTrackAtPad, Float_t yTrackAtPad, Int_t chamber)
 204 {
 205     //
 206     /// Decides if the digit belongs to the trigger track.
 207     //
 208
 209     Float_t numOfHalfWidth = 5.;
 210     Bool_t match = kFALSE;
 211     Float_t maxDistX = dpx;
 212     if(fReproduceTrigResponse && chamber>=2) maxDistX = 3.*dpx;// Non-bending plane: check the +- 1 strip between stations
 213     if(!fReproduceTrigResponse)maxDistX = numOfHalfWidth*dpx;
 214     Float_t maxDistY = dpy;
 215     if(fReproduceTrigResponse && chamber%2) maxDistY = 3*dpy;// bending plane: check the +- 1 strip between planes in the same station
 216     if(!fReproduceTrigResponse) maxDistY = numOfHalfWidth*dpy;
 217     Float_t deltaX = TMath::Abs(xPad-xTrackAtPad);
 218     Float_t deltaY = TMath::Abs(yPad-yTrackAtPad);
 219     if(deltaX<=maxDistX && deltaY<=maxDistY)match = kTRUE;
 220     return match;
 221 }
 222
 223
 224 //_____________________________________________________________________________
 225 Bool_t AliMUONTriggerChamberEff::IsDiffLocalBoard(Int_t currDetElemId, Int_t iy, Int_t detElemIdP1, Int_t iyDigitP1) const
 226 {
 227     //
 228     /// Determins if the digits belong to the same local board.
 229     /// Used only if one wants to reproduce the trigger algorithm result.
 230     /// (fReproduceTrigResponse = kTRUE).
 231     //
 232
 233     Bool_t isDiff = kTRUE;
 234     if(detElemIdP1<0 || iyDigitP1<0)return kFALSE;
 235     Int_t currSlat = currDetElemId%100;
 236     Int_t slatP1 = detElemIdP1%100;
 237     Int_t currLoc = iy/16;
 238     Int_t locP1 = iyDigitP1/16;
 239     if(currSlat==slatP1 && currLoc==locP1)isDiff = kFALSE;
 240     return isDiff;
 241 }
 242
 243
 244 //_____________________________________________________________________________
 245 void AliMUONTriggerChamberEff::CalculateEfficiency(Int_t trigger44, Int_t trigger34,
 246                                                    Float_t &efficiency, Float_t &error, Bool_t failuresAsInput)
 247 {
 248     //
 249     /// Returns the efficiency.
 250     //
 251
 252     efficiency=-9.;
 253     error=0.;
 254     if(trigger34>0){
 255         efficiency=(Double_t)trigger44/((Double_t)trigger34);
 256         if(failuresAsInput)efficiency=1.-(Double_t)trigger44/((Double_t)trigger34);
 257     }
 258     Double_t q = TMath::Abs(1-efficiency);
 259     if(efficiency<0)error=0.0;
 260     else error = TMath::Sqrt(efficiency*q/((Double_t)trigger34));
 261 }
 262
 263
 264 //_____________________________________________________________________________
 265 Int_t AliMUONTriggerChamberEff::DetElemIdFromPos(Float_t x, Float_t y, Int_t chamber, Int_t cathode)
 266 {
 267     //
 268     /// Given the (x,y) position in the chamber,
 269     /// it returns the corresponding slat
 270     //
 271
 272     Int_t resultingDetElemId = -1;
 273     AliMpDEIterator it;
 274     const AliMUONGeometryTransformer *kGeomTransformer = fMUON->GetGeometryTransformer();
 275     AliMUONSegmentation *segmentation = fMUON->GetSegmentation();
 276     for ( it.First(chamber-1); ! it.IsDone(); it.Next() ){
 277         Int_t detElemId = it.CurrentDEId();
 278
 279         if (  segmentation->HasDE(detElemId) ){
 280             const AliMpVSegmentation* seg =
 281                 AliMpSegmentation::Instance()
 282                   ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 283             if (seg){
 284                 Float_t deltax = seg->Dimensions().X();
 285                 Float_t deltay = seg->Dimensions().Y();
 286                 Float_t xlocal1 =  -deltax;
 287                 Float_t ylocal1 =  -deltay;
 288                 Float_t xlocal2 =  +deltax;
 289                 Float_t ylocal2 =  +deltay;
 290                 Float_t xg01, yg01, zg1, xg02, yg02, zg2;
 291                 kGeomTransformer->Local2Global(detElemId, xlocal1, ylocal1, 0, xg01, yg01, zg1);
 292                 kGeomTransformer->Local2Global(detElemId, xlocal2, ylocal2, 0, xg02, yg02, zg2);
 293
 294                 Float_t xg1 = xg01, xg2 = xg02, yg1 = yg01, yg2 = yg02;
 295
 296                 if(xg01>xg02){
 297                     xg1 = xg02;
 298                     xg2 = xg01;
 299                 }
 300                 if(yg01>yg02){
 301                     yg1 = yg02;
 302                     yg2 = yg01;
 303                 }
 304
 305                 if(x>=xg1 && x<=xg2 && y>=yg1 && y<=yg2){
 306                     resultingDetElemId = detElemId;
 307                     break;
 308                 }
 309             }
 310         }
 311     }
 312     return resultingDetElemId;
 313 }
 314
 315
 316 //_____________________________________________________________________________
 317 void AliMUONTriggerChamberEff::LocalBoardFromPos(Float_t x, Float_t y, Int_t detElemId, Int_t cathode, Int_t localBoard[4])
 318 {
 319     //
 320     /// Given the (x,y) position in the chamber,
 321     /// it returns the corresponding local board
 322     //
 323
 324     for(Int_t loc=0; loc<4; loc++){
 325         localBoard[loc]=-1;
 326     }
 327     const AliMUONGeometryTransformer *kGeomTransformer = fMUON->GetGeometryTransformer();
 328     Float_t xl, yl, zl;
 329     kGeomTransformer->Global2Local(detElemId, x, y, 0, xl, yl, zl);
 330     TVector2 pos(xl,yl);
 331     const AliMpVSegmentation* seg =
 332         AliMpSegmentation::Instance()
 333           ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 334     if (seg){
 335         AliMpPad pad = seg->PadByPosition(pos,kFALSE);
 336         for (Int_t loc=0; loc<pad.GetNofLocations(); loc++){
 337             AliMpIntPair location = pad.GetLocation(loc);
 338             localBoard[loc] = location.GetFirst();
 339         }
 340     }
 341 }
 342
 343
 344 //_____________________________________________________________________________
 345 void AliMUONTriggerChamberEff::PrintTrigger(AliMUONGlobalTrigger *globalTrig)
 346 {
 347     //
 348     /// Print trigger response.
 349     //
 350
 351     printf("===================================================\n");
 352     printf(" Global Trigger output\t \tLow pt\tHigh pt\n");
 353
 354     printf(" number of Single:\t \t");
 355     printf("%i\t",globalTrig->SingleLpt());
 356     printf("%i\t",globalTrig->SingleHpt());
 357     printf("\n");
 358
 359     printf(" number of UnlikeSign pair:\t");
 360     printf("%i\t",globalTrig->PairUnlikeLpt());
 361     printf("%i\t",globalTrig->PairUnlikeHpt());
 362     printf("\n");
 363
 364     printf(" number of LikeSign pair:\t");
 365     printf("%i\t",globalTrig->PairLikeLpt());
 366     printf("%i\t",globalTrig->PairLikeHpt());
 367     printf("\n");
 368     printf("===================================================\n");
 369     printf("\n");
 370 }
 371
 372 void AliMUONTriggerChamberEff::PerformTriggerChamberEff(const char* outputDir)
 373 {
 374     //
 375     /// Main method.
 376     /// It loops over the the trigger rec. tracks in the event.
 377     /// Then it search for matching digits around the track.
 378     /// Finally it calculates the efficiency for each trigger board.
 379     /// Files with calculated efficiency are placed in the user defined outputDir.
 380     //
 381
 382     enum {kBending, kNonBending};
 383     Int_t evtBeforePrint = 1000;
 384     Float_t rad2deg = 180./TMath::Pi();
 385
 386     Int_t chOrder[] = {0,2,1,3};
 387     Int_t station[] = {0,0,1,1};
 388     Float_t zRealMatch[fgkNchambers] = {0.0};
 389     Float_t correctFactor[fgkNcathodes] = {1.};
 390
 391     Bool_t match[fgkNchambers][fgkNcathodes] = {{kFALSE}};
 392
 393     TClonesArray *recTrigTracksArray = 0x0;
 394     AliMUONTriggerTrack *recTrigTrack = 0x0;
 395     AliMUONDigit * mDigit = 0x0;
 396
 397     Float_t zMeanChamber[fgkNchambers];
 398     for(Int_t ch=0; ch<fgkNchambers; ch++){
 399         zMeanChamber[ch] = AliMUONConstants::DefaultChamberZ(10+ch);
 400     }
 401
 402     TClonesArray * globalTrigger = 0x0;
 403     AliMUONGlobalTrigger * gloTrg = 0x0;
 404
 405     Int_t partNumOfTrig[fgkNchambers][fgkNcathodes] = {{0}};
 406     Int_t totNumOfTrig[fgkNchambers][fgkNcathodes] = {{0}};
 407     Int_t atLeast1MuPerEv[fgkNchambers][fgkNcathodes] = {{0}};
 408     Int_t digitPerTrack[fgkNcathodes] = {0};
 409
 410     Float_t trackIntersectCh[2][fgkNchambers]={{0.0}};
 411
 412     Int_t slatInPlane1[2][fgkNcathodes];
 413     Int_t iyDigitInPlane1[2][fgkNcathodes];
 414
 415     const Int_t kMaxNumOfTracks = 10;
 416     Int_t trigScheme[kMaxNumOfTracks][fgkNchambers][fgkNcathodes]={{{0}}};
 417     Int_t triggeredDigits[kMaxNumOfTracks][fgkNchambers][fgkNcathodes] = {{{-1}}};
 418     Int_t slatThatTriggered[kMaxNumOfTracks][fgkNchambers][fgkNcathodes]={{{-1}}};
 419     Int_t boardThatTriggered[kMaxNumOfTracks][fgkNchambers][fgkNcathodes][4]={{{{-1}}}};
 420     Int_t nboard[4]={-1};
 421     Int_t ineffBoard[4]={-1};
 422
 423     const Int_t kMaxNumOfDigits = 20;
 424     Int_t detElOfDigitsInData[kMaxNumOfDigits][fgkNchambers][fgkNcathodes] = {{{-1}}};
 425
 426     char filename[150];
 427     FileStat_t fs;
 428
 429     if(fFirstRun<0)fFirstRun=fLastRun=-1;
 430
 431     for(Int_t iRun = fFirstRun; iRun <= fLastRun; iRun++){// Loop over runs
 432     // open run loader and load gAlice
 433     if(fFirstRun>=0){
 434         cout<<"\n\nRun = "<<iRun<<endl;
 435         sprintf(filename, "%s/run%i/galice.root", fGaliceDir.Data(), iRun);
 436         if(gSystem->GetPathInfo(filename,fs)){
 437             cout<<"Warning: "<<filename<<" not found. Skip to next one"<<endl;
 438             continue;
 439         }
 440         cout<<"Opening file "<<filename<<endl;
 441         SetGaliceFile(filename);
 442     }
 443
 444     for (Int_t ievent=fFirstEvent; ievent<=fLastEvent; ievent++) { // event loop
 445         Bool_t isClearEvent = kTRUE;
 446
 447         for(Int_t ch=0; ch<fgkNchambers; ch++){
 448             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 449                 partNumOfTrig[ch][cath]=0;
 450                 match[ch][cath]=kFALSE;
 451                 for(Int_t itrack=0; itrack<kMaxNumOfTracks; itrack++){
 452                     triggeredDigits[itrack][ch][cath]=-1;
 453                     slatThatTriggered[itrack][ch][cath]=-1;
 454                     for(Int_t loc=0; loc<4; loc++){
 455                         boardThatTriggered[itrack][ch][cath][loc]=-1;
 456                     }
 457                 }
 458                 for(Int_t idig=0; idig<kMaxNumOfDigits; idig++){
 459                     detElOfDigitsInData[idig][ch][cath]=-1;
 460                 }
 461             }
 462         }
 463
 464         fRunLoader->GetEvent(ievent);
 465         if (ievent%evtBeforePrint==0) printf("\t Event = %d\n",ievent);
 466
 467         fData->SetTreeAddress("RL");
 468         fData->GetRecTriggerTracks();
 469         recTrigTracksArray = fData->RecTriggerTracks();
 470         Int_t nRecTrigTracks = (Int_t) recTrigTracksArray->GetEntriesFast();
 471
 472         fData->SetTreeAddress("D,GLT");
 473         fData->GetDigits();
 474
 475         const AliMUONGeometryTransformer* kGeomTransformer = fMUON->GetGeometryTransformer();
 476
 477         for (Int_t iRecTrigTrack=0; iRecTrigTrack<nRecTrigTracks; iRecTrigTrack++) {
 478             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 479                 digitPerTrack[cath]=0;
 480                 for(Int_t sta=0; sta<2; sta++){
 481                     slatInPlane1[sta][cath] = -9999;
 482                     iyDigitInPlane1[sta][cath] = -9999;
 483                 }
 484             }
 485
 486             Bool_t doubleCountTrack = kFALSE;
 487
 488             // reading info from tracks
 489             recTrigTrack = (AliMUONTriggerTrack *)recTrigTracksArray->At(iRecTrigTrack);
 490             Float_t x11 = recTrigTrack->GetX11();// x position (info from non-bending plane)
 491             Float_t y11 = recTrigTrack->GetY11();// y position (info from bending plane)
 492             Float_t thetaX = recTrigTrack->GetThetax();
 493             Float_t thetaY = recTrigTrack->GetThetay();
 494
 495             if(fDebugLevel>=3)printf("\tEvent = %i, Track = %i\npos from track: (x,y) = (%f, %f), (thetaX, thetaY) = (%f, %f)\n",ievent,iRecTrigTrack,x11,y11,thetaX*rad2deg,thetaY*rad2deg);
 496
 497             for(Int_t ch=0; ch<fgkNchambers; ch++) {
 498                 zRealMatch[ch] = zMeanChamber[ch];
 499                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 500                     trigScheme[iRecTrigTrack][ch][cath] = 0;
 501                 }
 502             }
 503
 504             for(Int_t ch=0; ch<fgkNchambers; ch++) { // chamber loop
 505                 Int_t currCh = chOrder[ch];
 506                 Int_t ichamber = 10+currCh;
 507                 Int_t currStation = station[currCh];
 508                 TClonesArray* digits = fData->Digits(ichamber);
 509                 digits->Sort();
 510                 Int_t ndigits = (Int_t)digits->GetEntriesFast();
 511                 if(fDebugLevel>=2){
 512                     if(fDebugLevel<3)printf("\tEvent = %i, Track = %i\n", ievent, iRecTrigTrack);
 513                     printf("DigitNum: %i digits detected\n",ndigits);
 514                 }
 515
 516                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 517                     correctFactor[cath]=1.;
 518                 }
 519                 // calculate corrections to trigger track theta
 520                 if(ch>=1)correctFactor[kNonBending] = zMeanChamber[0]/zRealMatch[0];// corrects x position
 521                 if(ch>=2)correctFactor[kBending] = (zMeanChamber[2] - zMeanChamber[0]) / (zRealMatch[2] - zRealMatch[0]);// corrects y position
 522
 523                 // searching track intersection with chambers (first approximation)
 524                 Float_t deltaZ = zMeanChamber[currCh] - zMeanChamber[0];
 525                 trackIntersectCh[0][currCh] = zMeanChamber[currCh] * TMath::Tan(thetaX) * correctFactor[kNonBending];// x position (info from non-bending plane)
 526                 trackIntersectCh[1][currCh] = y11 + deltaZ * TMath::Tan(thetaY) * correctFactor[kBending];// y position (info from bending plane)
 527
 528                 for(Int_t idigit=0; idigit<ndigits; idigit++) { // digit loop
 529                     mDigit = (AliMUONDigit*)digits->At(idigit);
 530                     for(Int_t loc=0; loc<4; loc++){
 531                         nboard[loc]=-1;
 532                     }
 533
 534                     // searching loaded digit global position and dimension
 535                     Int_t detElemId = mDigit->DetElemId();
 536                     Int_t cathode = mDigit->Cathode();
 537                     Int_t ix = mDigit->PadX();
 538                     Int_t iy = mDigit->PadY();
 539                     Float_t xpad, ypad, zpad;
 540                     if(detElOfDigitsInData[idigit][ch][cathode]==-1)detElOfDigitsInData[idigit][ch][cathode] = detElemId;
 541
 542                     if(fDebugLevel>=2)printf("cathode = %i\n",cathode);
 543                     const AliMpVSegmentation* seg = AliMpSegmentation::Instance()
 544                       ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 545
 546                     AliMpPad pad = seg->PadByIndices(AliMpIntPair(ix,iy),kTRUE);
 547                     for (Int_t loc=0; loc<pad.GetNofLocations(); loc++){
 548                         AliMpIntPair location = pad.GetLocation(loc);
 549                         nboard[loc] = location.GetFirst();
 550                     }
 551
 552                     // get the pad position and dimensions
 553                     Float_t xlocal1 = pad.Position().X();
 554                     Float_t ylocal1 = pad.Position().Y();
 555                     Float_t dpx = pad.Dimensions().X();
 556                     Float_t dpy = pad.Dimensions().Y();
 557
 558                     kGeomTransformer->Local2Global(detElemId, xlocal1, ylocal1, 0, xpad, ypad, zpad);
 559
 560                     if(fDebugLevel>=3)printf("ch = %i\t cath = %i\tpad = (%4.1f, %4.1f, %4.1f)\tsize = (%3.1f, %3.1f)\n",currCh,cathode,xpad,ypad,zpad,dpx,dpy);
 561
 562                     // searching track intersection with chambers (second approximation)
 563                     if(ch>=2){
 564                         deltaZ = zpad - zRealMatch[0];
 565                         trackIntersectCh[1][currCh] = y11 + deltaZ * TMath::Tan(thetaY) * correctFactor[kBending];// y position (info from bending plane)
 566                     }
 567
 568                     // deciding if digit matches track
 569                     Bool_t isDiffLocBoard = kFALSE;
 570                     if(fReproduceTrigResponse)isDiffLocBoard = IsDiffLocalBoard(detElemId, iy, slatInPlane1[currStation][cathode], iyDigitInPlane1[currStation][cathode]);
 571                     Bool_t matchPad = PadMatchTrack(xpad, ypad, dpx, dpy, trackIntersectCh[0][currCh], trackIntersectCh[1][currCh], currCh);
 572
 573                     if(matchPad && ch<2){
 574                         slatInPlane1[currStation][cathode] = detElemId;
 575                         iyDigitInPlane1[currStation][cathode] = iy;
 576                         if(fDebugLevel>=3)printf("slatInPlane1[%i][%i] = %i\tiyDigitInPlane1[%i][%i] = %i\n",currStation,cathode,slatInPlane1[currStation][cathode],currStation,cathode,iyDigitInPlane1[currStation][cathode]);
 577                     }
 578
 579                     if(isDiffLocBoard && fDebugLevel>=1)printf("\tDifferent local board\n");
 580
 581                     match[currCh][cathode] = (matchPad && !isDiffLocBoard);
 582
 583                     if(match[currCh][cathode]){
 584                         digitPerTrack[cathode]++;
 585                         trigScheme[iRecTrigTrack][currCh][cathode]++;
 586                         triggeredDigits[iRecTrigTrack][currCh][cathode] = idigit;
 587                         slatThatTriggered[iRecTrigTrack][currCh][cathode] = detElemId;
 588                         for(Int_t loc=0; loc<4; loc++){
 589                             boardThatTriggered[iRecTrigTrack][currCh][cathode][loc] = nboard[loc];
 590                         }
 591                         if(digitPerTrack[cathode]>4 && !fReproduceTrigResponse)isClearEvent = kFALSE;
 592                     }
 593
 594                     // in case of match, store real z position of the chamber
 595                     if(matchPad)zRealMatch[currCh] = zpad;
 596
 597                 } // end digit loop
 598             } // end chamber loop
 599
 600             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 601                 if(digitPerTrack[cath]<3 && !fReproduceTrigResponse)isClearEvent = kFALSE;
 602             }
 603
 604             if(!isClearEvent && !fReproduceTrigResponse){
 605                 fData->ResetDigits();
 606                 continue;
 607             }
 608
 609             Int_t commonDigits = 0;
 610             Int_t doubleTrack = -1;
 611             for(Int_t itrack=0; itrack<iRecTrigTrack; itrack++){
 612                 for(Int_t ch=0; ch<fgkNchambers; ch++){
 613                     if(triggeredDigits[itrack][ch][kBending]==triggeredDigits[iRecTrigTrack][ch][kBending])commonDigits++;
 614                 }
 615                 if(commonDigits>=2){
 616                     doubleCountTrack=kTRUE;
 617                     doubleTrack = itrack;
 618                     break;
 619                 }
 620             }
 621
 622             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 623                 if(!doubleCountTrack || fReproduceTrigResponse){
 624                     Int_t is44 = 1;
 625                     Bool_t goodForSlatEff=kTRUE;
 626                     Bool_t goodForBoardEff=kTRUE;
 627                     Int_t firstSlat=slatThatTriggered[iRecTrigTrack][0][cath]%100;
 628                     if(firstSlat<0)firstSlat=slatThatTriggered[iRecTrigTrack][1][cath]%100;
 629                     Int_t firstBoard=boardThatTriggered[iRecTrigTrack][0][0][0];
 630                     if(firstBoard<0)firstBoard=boardThatTriggered[iRecTrigTrack][1][0][0];
 631                     for(Int_t ch=0; ch<fgkNchambers; ch++){
 632                         is44 *= trigScheme[iRecTrigTrack][ch][cath];
 633                         Int_t currSlat=slatThatTriggered[iRecTrigTrack][ch][cath]%100;
 634                         if(currSlat<0)continue;
 635                         if(currSlat!=firstSlat)goodForSlatEff=kFALSE;
 636                         Bool_t atLeastOneLoc=kFALSE;
 637                         for(Int_t loc=0; loc<4; loc++){
 638                             Int_t currBoard = boardThatTriggered[iRecTrigTrack][ch][cath][loc];
 639                             if(currBoard==firstBoard){
 640                                 atLeastOneLoc=kTRUE;
 641                                 break;
 642                             }
 643                         }
 644                         if(!atLeastOneLoc)goodForBoardEff=kFALSE;
 645                     }
 646                     if(fDebugLevel==1)printf("\tEvent = %i, Track = %i\n", ievent, iRecTrigTrack);
 647                     if(is44==1){
 648                         fTrigger44[cath]++;
 649                         if(fDebugLevel>=1)printf("Trigger44[%i] = %i\n",cath,fTrigger44[cath]);
 650                         if(goodForSlatEff){
 651                             for(Int_t ch=0; ch<fgkNchambers; ch++){
 652                                 for(Int_t slat=0; slat<fgkNslats; slat++){
 653                                     Int_t corrDetEl = (ch+11)*100 + slat;
 654                                     if(corrDetEl==slatThatTriggered[iRecTrigTrack][ch][cath]){
 655                                         fHitPerSlat[ch][cath][slat]++;
 656                                         if(fDebugLevel>=1)printf("Slat that triggered = %i\n",corrDetEl);
 657                                         if(goodForBoardEff && firstBoard>0){
 658                                             fHitPerBoard[ch][cath][firstBoard-1]++;
 659                                             if(fDebugLevel>=1)printf("Board that triggered = %i\n",firstBoard);
 660                                         }
 661                                         else if(fDebugLevel>=1)printf("Event = %i, Track = %i: Particle crossed different boards: rejected!\n",ievent,iRecTrigTrack);
 662                                     }
 663                                 }
 664                             }
 665                         }
 666                         else printf("Event = %i, Track = %i: Particle crossed different slats: rejected!\n",ievent,iRecTrigTrack);
 667                     }
 668                     if(digitPerTrack[cath]==3){
 669                         for(Int_t ch=0; ch<fgkNchambers; ch++){
 670                             if(match[ch][cath])partNumOfTrig[ch][cath]++;
 671                             if(trigScheme[iRecTrigTrack][ch][cath]==0){
 672                                 fTrigger34[ch][cath]++;
 673                                 if(fDebugLevel>=1)printf("Trigger34[%i][%i] = %i\n",ch,cath,fTrigger34[ch][cath]);
 674                                 if(!goodForSlatEff){
 675                                     printf("Event %i, Track = %i: Particle crossed different slats: rejected!\n",ievent,iRecTrigTrack);
 676                                     continue;
 677                                 }
 678                                 Int_t ineffSlat = DetElemIdFromPos(trackIntersectCh[0][ch], trackIntersectCh[1][ch], 11+ch, cath);
 679                                 if(fDebugLevel>=1)printf("Slat non efficient = %i\n",ineffSlat);
 680                                 if(ineffSlat>0){
 681                                     Int_t slatInCh = ineffSlat%100;
 682                                     fInefficientSlat[ch][cath][slatInCh]++;
 683                                     for(Int_t idig=0; idig<kMaxNumOfDigits; idig++){
 684                                         if(ineffSlat==detElOfDigitsInData[idig][ch][cath])cout<<"Warning: "<<ineffSlat<<" is not inefficient!!!"<<endl;
 685                                     }
 686                                     LocalBoardFromPos(trackIntersectCh[0][ch], trackIntersectCh[1][ch], ineffSlat, cath, ineffBoard);
 687                                     Int_t boardNonEff=-1;
 688                                     for(Int_t loc=0; loc<4; loc++){
 689                                         if(ineffBoard[loc]==firstBoard){
 690                                             boardNonEff=ineffBoard[loc];
 691                                             break;
 692                                         }
 693                                     }
 694                                     if(fDebugLevel>=1)printf("Board non efficient = %i\n",boardNonEff);
 695                                     if(boardNonEff>0)fInefficientBoard[ch][cath][boardNonEff-1]++;
 696                                     else if(fDebugLevel>=1){
 697                                         printf("Inefficient board should be %i.\tBoards found:\n", firstBoard);
 698                                         for(Int_t loc=0; loc<4; loc++){
 699                                             printf("%i\t",ineffBoard[loc]);
 700                                         }
 701                                         printf("\n");
 702                                     }
 703
 704                                 }
 705                             }
 706                         }
 707                     }
 708                 }
 709                 else if(doubleCountTrack){
 710                     if(fDebugLevel<=1)printf("\n\tEvent = %i, Track = %i: ", ievent,iRecTrigTrack);
 711                     printf("Double Count Track: %i similar to %i. Track rejected!\n",iRecTrigTrack, doubleTrack);
 712                 }
 713             }
 714         }// end trigger tracks loop
 715         if(nRecTrigTracks<=0){
 716             fData->ResetDigits();
 717             continue;
 718         }
 719         for(Int_t ch=0; ch<fgkNchambers; ch++){
 720             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 721                 totNumOfTrig[ch][cath] += partNumOfTrig[ch][cath];
 722                 if(partNumOfTrig[ch][cath]>0)atLeast1MuPerEv[ch][cath]++;
 723             }
 724         }
 725
 726         if(fPrintInfo){
 727             //Global trigger
 728             globalTrigger = fData->GlobalTrigger();
 729             Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
 730
 731             for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
 732                 gloTrg = (AliMUONGlobalTrigger*)globalTrigger->At(iglobal);
 733             }
 734             PrintTrigger(gloTrg);
 735             InfoDigit();
 736             cout<<"\n"<<endl;
 737         }
 738
 739         fData->ResetDigits();
 740     }// end event loop
 741     if(fFirstRun>=0)CleanGalice();
 742     } //end loop over run
 743
 744     // Write output data
 745     WriteEfficiencyMap(outputDir);
 746
 747     WriteOutput(outputDir, totNumOfTrig, atLeast1MuPerEv);
 748 }
 749
 750 //_____________________________________________________________________________
 751 void AliMUONTriggerChamberEff::WriteOutput(const char* outputDir, Int_t totNumOfTrig[4][2], Int_t atLeast1MuPerEv[4][2])
 752 {
 753     //
 754     /// Writes information on calculated efficiency.
 755     /// It writes: triggerChamberEff.root file containing efficiency histograms.
 756     ///
 757     /// In addition a text file triggerChamberEff.out is created,
 758     /// with further informations on efficiencies.
 759     //
 760
 761     char *cathodeName[fgkNcathodes]={"Bending plane", "Non-Bending plane"};
 762     char *cathCode[fgkNcathodes] = {"bendPlane", "nonBendPlane"};
 763
 764     char outFileName[100];
 765     sprintf(outFileName, "%s/triggerChamberEff.out",outputDir);
 766     FILE *outfile = fopen(outFileName, "w");
 767     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 768         fprintf(outfile,"%s:\n",cathodeName[cath]);
 769         for(Int_t ch=0; ch<fgkNchambers; ch++){
 770             fprintf(outfile,"Total number of muon triggers chamber 1%i = %i\n",ch+1,totNumOfTrig[ch][cath]);
 771         }
 772         fprintf(outfile,"\n");
 773     }
 774     fprintf(outfile,"\n");
 775     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 776         fprintf(outfile,"%s:\n",cathodeName[cath]);
 777         for(Int_t ch=0; ch<fgkNchambers; ch++){
 778             fprintf(outfile,"At least 1 muon triggered chamber 1%i = %i\n",ch+1, atLeast1MuPerEv[ch][cath]);
 779         }
 780         fprintf(outfile,"\n");
 781     }
 782     fprintf(outfile,"\n\n");
 783     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 784         fprintf(outfile,"%s:\n",cathodeName[cath]);
 785         fprintf(outfile,"Number of triggers where all chambers counted = %i\n",fTrigger44[cath]);
 786         fprintf(outfile,"\n");
 787     }
 788     fprintf(outfile,"\n");
 789     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 790         fprintf(outfile,"%s:\n",cathodeName[cath]);
 791         for(Int_t ch=0; ch<fgkNchambers; ch++){
 792             fprintf(outfile,"Number of triggers where chamber 1%i did not count = %i\n",ch+1,fTrigger34[ch][cath]);
 793         }
 794         fprintf(outfile,"\n");
 795     }
 796
 797     fprintf(outfile,"\n");
 798     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 799         fprintf(outfile,"%s:\n",cathodeName[cath]);
 800         for(Int_t ch=0; ch<fgkNchambers; ch++){
 801             Int_t sumIneff = 0, sumHits = 0;
 802             fprintf(outfile,"\n Chamber %1i\n", ch+1);
 803             for(Int_t slat=0; slat<fgkNslats; slat++){
 804                 fprintf(outfile,"Number of triggers where slat %2i - did not count = %5i - was hit (hit%sCh%iSlat%i) = %5i\n",slat,fInefficientSlat[ch][cath][slat],cathCode[cath],11+ch,slat,fHitPerSlat[ch][cath][slat]);
 805                 sumIneff += fInefficientSlat[ch][cath][slat];
 806                 sumHits += fHitPerSlat[ch][cath][slat];
 807             }
 808             fprintf(outfile,"Number of triggers where chamber %1i - did not count = %5i - was hit (hit%sCh%i) = %5i\n",ch+1,sumIneff,cathCode[cath],11+ch,sumHits);
 809         }
 810         fprintf(outfile,"\n");
 811     }
 812     fclose(outfile);
 813
 814     sprintf(outFileName, "%s/triggerChamberEff.root",outputDir);
 815     TFile *outputHistoFile = new TFile(outFileName,"RECREATE");
 816     TDirectory *dir = gDirectory;
 817
 818     enum {kSlatIn11, kSlatIn12, kSlatIn13, kSlatIn14, kChamberEff};
 819     char *yAxisTitle = "trigger efficiency (a.u.)";
 820     char *xAxisTitle = "chamber";
 821
 822     TH1F *histo[fgkNcathodes][fgkNchambers+1];
 823     TH1F *histoBoard[fgkNcathodes][fgkNchambers];
 824
 825     char histoName[30];
 826     char histoTitle[90];
 827
 828     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 829         for(Int_t ch=0; ch<fgkNchambers+1; ch++){
 830             if(ch==kChamberEff){
 831                 sprintf(histoName, "%sChamberEff", cathCode[cath]);
 832                 sprintf(histoTitle, "Chamber efficiency %s", cathCode[cath]);
 833                 histo[cath][ch] = new TH1F(histoName, histoTitle, fgkNchambers, 11-0.5, 15-0.5);
 834                 histo[cath][ch]->SetXTitle(xAxisTitle);
 835                 histo[cath][ch]->SetYTitle(yAxisTitle);
 836                 histo[cath][ch]->GetXaxis()->SetNdivisions(fgkNchambers);
 837             }
 838             else {
 839                 sprintf(histoName, "%sSlatEffChamber%i", cathCode[cath], 11+ch);
 840                 sprintf(histoTitle, "Chamber %i: slat efficiency %s", 11+ch, cathCode[cath]);
 841                 histo[cath][ch] = new TH1F(histoName, histoTitle, fgkNslats, 0-0.5, fgkNslats-0.5);
 842                 histo[cath][ch]->SetXTitle("slat");
 843                 histo[cath][ch]->SetYTitle(yAxisTitle);
 844                 histo[cath][ch]->GetXaxis()->SetNdivisions(fgkNslats);
 845
 846                 sprintf(histoName, "%sBoardEffChamber%i", cathCode[cath], 11+ch);
 847                 sprintf(histoTitle, "Chamber %i: board efficiency %s", 11+ch, cathCode[cath]);
 848                 histoBoard[cath][ch] = new TH1F(histoName, histoTitle, fgkNboards, 1-0.5, fgkNboards+1.-0.5);
 849                 histoBoard[cath][ch]->SetXTitle("boards");
 850                 histoBoard[cath][ch]->SetYTitle(yAxisTitle);
 851                 histoBoard[cath][ch]->GetXaxis()->SetNdivisions(fgkNboards);
 852             }
 853         }
 854     }
 855
 856     Float_t efficiency, efficiencyError;
 857
 858     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 859         for(Int_t ch=0; ch<fgkNchambers; ch++){
 860             for(Int_t slat=0; slat<fgkNslats; slat++){
 861                 CalculateEfficiency(fHitPerSlat[ch][cath][slat], fHitPerSlat[ch][cath][slat]+fInefficientSlat[ch][cath][slat], efficiency, efficiencyError, kFALSE);
 862                 histo[cath][ch]->SetBinContent(slat+1, efficiency);
 863                 histo[cath][ch]->SetBinError(slat+1, efficiencyError);
 864             }
 865             CalculateEfficiency(fTrigger44[cath], fTrigger34[ch][cath]+fTrigger44[cath], efficiency, efficiencyError, kFALSE);
 866             histo[cath][kChamberEff]->SetBinContent(ch+1, efficiency);
 867             histo[cath][kChamberEff]->SetBinError(ch+1, efficiencyError);
 868
 869             for(Int_t board=0; board<fgkNboards; board++){
 870                 CalculateEfficiency(fHitPerBoard[ch][cath][board], fHitPerBoard[ch][cath][board]+fInefficientBoard[ch][cath][board], efficiency, efficiencyError, kFALSE);
 871                 histoBoard[cath][ch]->SetBinContent(board+1, efficiency);
 872                 histoBoard[cath][ch]->SetBinError(board+1, efficiencyError);
 873             }
 874         }
 875     }
 876
 877 // write all histos
 878     outputHistoFile->cd();
 879     dir->GetList()->Write();
 880     outputHistoFile->Close();
 881 }
 882
 883
 884 //_____________________________________________________________________________
 885 void AliMUONTriggerChamberEff::WriteEfficiencyMap(const char* outputDir)
 886 {
 887     //
 888     /// Writes the calculated efficiency in the text file efficiencyCells.dat
 889     ///
 890     /// The file can be further put in $ALICE_ROOT/MUON/data
 891     /// and used to run simulations with measured trigger chamber efficiencies.
 892     //
 893
 894     Int_t effOutWidth=4;
 895
 896     Float_t efficiency, efficiencyError;
 897
 898     Int_t aCapo[] = {16, 38, 60, 76, 92, 108, 117, 133, 155, 177, 193, 209, 225, 234};
 899
 900     char filename[70];
 901     sprintf(filename, "%s/efficiencyCells.dat", outputDir);
 902     ofstream outFile(filename);
 903     outFile << "localBoards" << endl;
 904     for(Int_t ch=0; ch<fgkNchambers; ch++){
 905         //Print information
 906         outFile << "\n\ndetElemId:\t" << 11+ch;
 907         outFile << "00";
 908         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 909             outFile << "\n cathode:\t" << cath << endl;
 910             Int_t currLine=0;
 911             for(Int_t board=0; board<fgkNboards; board++){
 912
 913                 if(board==aCapo[currLine]){
 914                     outFile << endl;
 915                     currLine++;
 916                 }
 917                 CalculateEfficiency(fHitPerBoard[ch][cath][board], fHitPerBoard[ch][cath][board]+fInefficientBoard[ch][cath][board], efficiency, efficiencyError, kFALSE);
 918                 outFile << " " << setw(effOutWidth) << efficiency;
 919             }// loop on boards
 920             outFile << endl;
 921         }// loop on cathodes
 922     }// loop on chambers
 923 }
 924