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 "AliMUONRecData.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 }
  94
  95 //_____________________________________________________________________________
  96 AliMUONTriggerChamberEff::~AliMUONTriggerChamberEff()
  97 {
  98 /// Destructor
  99     fRunLoader->UnloadAll();
 100     delete fRunLoader;
 101     delete fData;
 102 }
 103
 104 //_____________________________________________________________________________
 105 void AliMUONTriggerChamberEff::SetGaliceFile(const char *galiceFile)
 106 {
 107     //
 108     /// Opens the galice.root and loads tracks and digits.
 109     //
 110
 111     fRunLoader = AliRunLoader::Open(galiceFile,"MUONFolder","READ");
 112     if (!fRunLoader)
 113     {
 114         AliError(Form("Error opening %s file \n",galiceFile));
 115     }
 116     else
 117     {
 118         fRunLoader->LoadgAlice();
 119         gAlice = fRunLoader->GetAliRun();
 120         fMUON = (AliMUON*)gAlice->GetModule("MUON");
 121
 122         if(fLastEvent<=0 || fLastEvent>fRunLoader->GetNumberOfEvents())fLastEvent = fRunLoader->GetNumberOfEvents()-1;
 123         if(fFirstEvent<0)fFirstEvent=0;
 124
 125
 126         AliLoader* loader = fRunLoader->GetLoader("MUONLoader");
 127         if ( loader )
 128         {
 129             fData = new AliMUONRecData(loader,"MUON","MUON");
 130             loader->LoadTracks("READ");
 131             loader->LoadDigits("READ");
 132         }
 133         else
 134         {
 135             AliError(Form("Could get MUONLoader"));
 136         }
 137     }
 138 }
 139
 140 //_____________________________________________________________________________
 141 void AliMUONTriggerChamberEff::CleanGalice()
 142 {
 143     //
 144     /// Unload all loaded data
 145     //
 146     delete fData;
 147     fData = NULL;
 148     fRunLoader->UnloadAll("all");
 149     delete fRunLoader;
 150     fRunLoader = NULL;
 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     fData->ResetTrigger();
 199 }
 200
 201
 202 //_____________________________________________________________________________
 203 Int_t AliMUONTriggerChamberEff::MatchingPad(Int_t &detElemId, Float_t coor[2], const AliMUONGeometryTransformer *kGeomTransformer, Bool_t isMatch[fgkNcathodes], Int_t nboard[fgkNcathodes][4], Float_t zRealMatch[fgkNchambers], Float_t y11)
 204 {
 205     //
 206     /// Check slat and board number of digit matching track
 207     //
 208
 209     enum {kBending, kNonBending};
 210
 211     Float_t minMatchDist[fgkNcathodes];
 212
 213     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 214         isMatch[cath]=kFALSE;
 215         minMatchDist[cath]=9999.;
 216     }
 217     Int_t iChamber = detElemId/100-1;
 218     Int_t ch = iChamber-10;
 219     Float_t oldDeltaZ = AliMUONConstants::DefaultChamberZ(iChamber) - AliMUONConstants::DefaultChamberZ(10);
 220     Float_t y = coor[1];
 221     Int_t iSlat = detElemId%100;
 222     Int_t trigDigitBendPlane = -1;
 223     TClonesArray* digits = fData->Digits(iChamber);
 224     digits->Sort();
 225     Int_t foundDetElemId = detElemId;
 226     Float_t foundZmatch=999.;
 227     Float_t yCoorAtPadZ=999.;
 228     Int_t ndigits = (Int_t)digits->GetEntriesFast();
 229     AliMUONDigit * mDigit = 0x0;
 230     for(Int_t idigit=0; idigit<ndigits; idigit++) { // digit loop
 231         mDigit = (AliMUONDigit*)digits->At(idigit);
 232         Int_t currDetElemId = mDigit->DetElemId();
 233         Int_t currSlat = currDetElemId%100;
 234         if(TMath::Abs(currSlat%18-iSlat%18)>1)continue; // Check neighbour slats
 235         Int_t cathode = mDigit->Cathode();
 236         Int_t ix = mDigit->PadX();
 237         Int_t iy = mDigit->PadY();
 238         Float_t xpad, ypad, zpad;
 239         const AliMpVSegmentation* seg = AliMpSegmentation::Instance()
 240             ->GetMpSegmentation(currDetElemId,AliMp::GetCathodType(cathode));
 241
 242         AliMpPad pad = seg->PadByIndices(AliMpIntPair(ix,iy),kTRUE);
 243         Float_t xlocal1 = pad.Position().X();
 244         Float_t ylocal1 = pad.Position().Y();
 245         Float_t dpx = pad.Dimensions().X();
 246         Float_t dpy = pad.Dimensions().Y();
 247         kGeomTransformer->Local2Global(currDetElemId, xlocal1, ylocal1, 0, xpad, ypad, zpad);
 248         if(fDebugLevel>2)printf("DetElemId = %i\tCathode = %i\t(x,y) Pad = (%i,%i) = (%.2f,%.2f)\tDim = (%.2f,%.2f)\tTrack = (%.2f,%.2f)\n",currDetElemId,cathode,ix,iy,xpad,ypad,dpx,dpy,coor[0],coor[1]);
 249         // searching track intersection with chambers (second approximation)
 250         if(ch%2==1){
 251             Float_t deltaZ = zpad - zRealMatch[0];
 252             y = (coor[1]-y11)*deltaZ/oldDeltaZ + y11;
 253             if(fDebugLevel>=3 && TMath::Abs(y-coor[1])>0.1)printf("oldDeltaZ = %7.2f   newDeltaZ = %7.2f\toldY = %7.2f   new y = %7.2f\n",oldDeltaZ,deltaZ,coor[1],y);
 254         }
 255         Float_t matchDist = PadMatchTrack(xpad, ypad, dpx, dpy, coor[0], y, ch);
 256         if(matchDist>minMatchDist[cathode])continue;
 257         isMatch[cathode] = kTRUE;
 258         minMatchDist[cathode] = matchDist;
 259         foundDetElemId = currDetElemId;
 260         foundZmatch=zpad;
 261         yCoorAtPadZ=y;
 262         if(cathode==kBending)trigDigitBendPlane = idigit;
 263         for (Int_t loc=0; loc<pad.GetNofLocations(); loc++){
 264             AliMpIntPair location = pad.GetLocation(loc);
 265             nboard[cathode][loc] = location.GetFirst();
 266         }
 267         for(Int_t loc=pad.GetNofLocations(); loc<4; loc++){
 268             nboard[cathode][loc]=-1;
 269         }
 270     }
 271     if(isMatch[kBending] || isMatch[kNonBending]){
 272         detElemId = foundDetElemId;
 273         zRealMatch[ch] = foundZmatch;
 274         coor[1] = yCoorAtPadZ;
 275         if(fDebugLevel>2){
 276             Int_t whichCathode=kBending;
 277             if(!isMatch[kBending])whichCathode=kNonBending;
 278         }
 279     }
 280     return trigDigitBendPlane;
 281 }
 282
 283 //_____________________________________________________________________________
 284 Float_t AliMUONTriggerChamberEff::PadMatchTrack(Float_t xPad, Float_t yPad, Float_t dpx, Float_t dpy,
 285                                                Float_t xTrackAtPad, Float_t yTrackAtPad, Int_t chamber)
 286 {
 287     //
 288     /// Decides if the digit belongs to the trigger track.
 289     //
 290
 291     Float_t maxDist = 3.;//cm
 292
 293     Float_t matchDist = 99999.;
 294
 295     Float_t deltaX = TMath::Abs(xPad-xTrackAtPad)-dpx;
 296     Float_t deltaY = TMath::Abs(yPad-yTrackAtPad)-dpy;
 297     Float_t maxDistX = maxDist;
 298     Float_t maxDistY = maxDist;
 299
 300     if(fReproduceTrigResponse){
 301         maxDistX = dpx;
 302         maxDistY = dpy;
 303         deltaX = TMath::Abs(xPad-xTrackAtPad);
 304         deltaY = TMath::Abs(yPad-yTrackAtPad);
 305         if(dpx<dpy && chamber>=2) maxDistX = 3.*dpx;// Non-bending plane: check the +- 1 strip between stations
 306         if(dpy<dpx && chamber%2) maxDistY = 3.*dpy;// bending plane: check the +- 1 strip between planes in the same station
 307     }
 308
 309     if(deltaX<=maxDistX && deltaY<=maxDistY)matchDist = deltaX*deltaX + deltaY*deltaY;
 310     return matchDist;
 311 }
 312
 313
 314 //_____________________________________________________________________________
 315 void AliMUONTriggerChamberEff::CalculateEfficiency(Int_t trigger44, Int_t trigger34,
 316                                                    Float_t &efficiency, Float_t &error, Bool_t failuresAsInput)
 317 {
 318     //
 319     /// Returns the efficiency.
 320     //
 321
 322     efficiency=-9.;
 323     error=0.;
 324     if(trigger34>0){
 325         efficiency=(Double_t)trigger44/((Double_t)trigger34);
 326         if(failuresAsInput)efficiency=1.-(Double_t)trigger44/((Double_t)trigger34);
 327     }
 328     Double_t q = TMath::Abs(1-efficiency);
 329     if(efficiency<0)error=0.0;
 330     else error = TMath::Sqrt(efficiency*q/((Double_t)trigger34));
 331 }
 332
 333
 334 //_____________________________________________________________________________
 335 Int_t AliMUONTriggerChamberEff::DetElemIdFromPos(Float_t x, Float_t y, Int_t chamber, Int_t cathode)
 336 {
 337     //
 338     /// Given the (x,y) position in the chamber,
 339     /// it returns the corresponding slat
 340     //
 341
 342     Int_t resultingDetElemId = -1;
 343     AliMpDEIterator it;
 344     const AliMUONGeometryTransformer *kGeomTransformer = fMUON->GetGeometryTransformer();
 345     AliMUONSegmentation *segmentation = fMUON->GetSegmentation();
 346     Float_t minDist = 999.;
 347     for ( it.First(chamber-1); ! it.IsDone(); it.Next() ){
 348         Int_t detElemId = it.CurrentDEId();
 349         Int_t ich = detElemId/100-10;
 350         Float_t tolerance=0.2*((Float_t)ich);
 351         Float_t currDist=9999.;
 352
 353         if (  segmentation->HasDE(detElemId) ){
 354             const AliMpVSegmentation* seg =
 355                 AliMpSegmentation::Instance()
 356                   ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 357             if (seg){
 358                 Float_t deltax = seg->Dimensions().X();
 359                 Float_t deltay = seg->Dimensions().Y();
 360                 Float_t xlocal1 =  -deltax;
 361                 Float_t ylocal1 =  -deltay;
 362                 Float_t xlocal2 =  +deltax;
 363                 Float_t ylocal2 =  +deltay;
 364                 Float_t xg01, yg01, zg1, xg02, yg02, zg2;
 365                 kGeomTransformer->Local2Global(detElemId, xlocal1, ylocal1, 0, xg01, yg01, zg1);
 366                 kGeomTransformer->Local2Global(detElemId, xlocal2, ylocal2, 0, xg02, yg02, zg2);
 367
 368                 Float_t xg1 = xg01, xg2 = xg02, yg1 = yg01, yg2 = yg02;
 369
 370                 if(xg01>xg02){
 371                     xg1 = xg02;
 372                     xg2 = xg01;
 373                 }
 374                 if(yg01>yg02){
 375                     yg1 = yg02;
 376                     yg2 = yg01;
 377                 }
 378
 379                 if(x>=xg1-tolerance && x<=xg2+tolerance && y>=yg1-tolerance && y<=yg2+tolerance){ // takes into account errors in extrapolation
 380                     if(y<yg1) currDist = yg1-y;
 381                     else if(y>yg2) currDist = y-yg2;
 382                     if(currDist<minDist) {
 383                         resultingDetElemId = detElemId;
 384                         minDist=currDist;
 385                         continue;
 386                     }
 387                     resultingDetElemId = detElemId;
 388                     break;
 389                 }
 390             }
 391         }
 392     }
 393     return resultingDetElemId;
 394 }
 395
 396
 397 //_____________________________________________________________________________
 398 void AliMUONTriggerChamberEff::LocalBoardFromPos(Float_t x, Float_t y, Int_t detElemId, Int_t cathode, Int_t localBoard[4])
 399 {
 400     //
 401     /// Given the (x,y) position in the chamber,
 402     /// it returns the corresponding local board
 403     //
 404
 405     for(Int_t loc=0; loc<4; loc++){
 406         localBoard[loc]=-1;
 407     }
 408     const AliMUONGeometryTransformer *kGeomTransformer = fMUON->GetGeometryTransformer();
 409     Float_t xl, yl, zl;
 410     kGeomTransformer->Global2Local(detElemId, x, y, 0, xl, yl, zl);
 411     TVector2 pos(xl,yl);
 412     const AliMpVSegmentation* seg =
 413         AliMpSegmentation::Instance()
 414           ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 415     if (seg){
 416         AliMpPad pad = seg->PadByPosition(pos,kFALSE);
 417         for (Int_t loc=0; loc<pad.GetNofLocations(); loc++){
 418             AliMpIntPair location = pad.GetLocation(loc);
 419             localBoard[loc] = location.GetFirst();
 420         }
 421     }
 422 }
 423
 424
 425 //_____________________________________________________________________________
 426 void AliMUONTriggerChamberEff::PrintTrigger(AliMUONGlobalTrigger *globalTrig)
 427 {
 428     //
 429     /// Print trigger response.
 430     //
 431
 432     printf("===================================================\n");
 433     printf(" Global Trigger output\t \tLow pt\tHigh pt\n");
 434
 435     printf(" number of Single:\t \t");
 436     printf("%i\t",globalTrig->SingleLpt());
 437     printf("%i\t",globalTrig->SingleHpt());
 438     printf("\n");
 439
 440     printf(" number of UnlikeSign pair:\t");
 441     printf("%i\t",globalTrig->PairUnlikeLpt());
 442     printf("%i\t",globalTrig->PairUnlikeHpt());
 443     printf("\n");
 444
 445     printf(" number of LikeSign pair:\t");
 446     printf("%i\t",globalTrig->PairLikeLpt());
 447     printf("%i\t",globalTrig->PairLikeHpt());
 448     printf("\n");
 449     printf("===================================================\n");
 450     printf("\n");
 451 }
 452
 453 void AliMUONTriggerChamberEff::PerformTriggerChamberEff(const char* outputDir)
 454 {
 455     //
 456     /// Main method.
 457     /// It loops over the the trigger rec. tracks in the event.
 458     /// Then it search for matching digits around the track.
 459     /// Finally it calculates the efficiency for each trigger board.
 460     /// Files with calculated efficiency are placed in the user defined outputDir.
 461     //
 462
 463     enum {kBending, kNonBending};
 464     Int_t evtBeforePrint = 1000;
 465     Float_t rad2deg = 180./TMath::Pi();
 466
 467     Int_t chOrder[fgkNchambers] = {0,2,1,3};
 468     Float_t zRealMatch[fgkNchambers] = {0.0};
 469     Float_t correctFactor[fgkNcathodes] = {1.};
 470
 471     Bool_t match[fgkNchambers][fgkNcathodes] = {{kFALSE}};
 472     Bool_t matchPad[fgkNcathodes]={kFALSE};
 473
 474     TClonesArray *recTrigTracksArray = 0x0;
 475     AliMUONTriggerTrack *recTrigTrack = 0x0;
 476
 477     Float_t zMeanChamber[fgkNchambers];
 478     for(Int_t ch=0; ch<fgkNchambers; ch++){
 479         zMeanChamber[ch] = AliMUONConstants::DefaultChamberZ(10+ch);
 480     }
 481
 482     TClonesArray * globalTrigger = 0x0;
 483     AliMUONGlobalTrigger * gloTrg = 0x0;
 484
 485     Int_t partNumOfTrig[fgkNchambers][fgkNcathodes] = {{0}};
 486     Int_t totNumOfTrig[fgkNchambers][fgkNcathodes] = {{0}};
 487     Int_t atLeast1MuPerEv[fgkNchambers][fgkNcathodes] = {{0}};
 488     Int_t digitPerTrack[fgkNcathodes] = {0};
 489
 490     Float_t trackIntersectCh[fgkNchambers][2]={{0.0}};
 491
 492     Int_t triggeredDigits[2][fgkNchambers] = {{-1}};
 493
 494     Int_t trigScheme[fgkNchambers][fgkNcathodes]={{0}};
 495     Int_t slatThatTriggered[fgkNchambers][fgkNcathodes]={{-1}};
 496     Int_t boardThatTriggered[fgkNchambers][fgkNcathodes][4]={{{-1}}};
 497     Int_t nboard[fgkNcathodes][4]={{-1}};
 498     Int_t ineffBoard[4]={-1};
 499
 500     char filename[150];
 501     FileStat_t fs;
 502
 503     if(fFirstRun<0)fFirstRun=fLastRun=-1;
 504
 505     for(Int_t iRun = fFirstRun; iRun <= fLastRun; iRun++){// Loop over runs
 506     // open run loader and load gAlice
 507     if(fFirstRun>=0){
 508         cout<<"\n\nRun = "<<iRun<<endl;
 509         sprintf(filename, "%s/run%i/galice.root", fGaliceDir.Data(), iRun);
 510         if(gSystem->GetPathInfo(filename,fs)){
 511             cout<<"Warning: "<<filename<<" not found. Skip to next one"<<endl;
 512             continue;
 513         }
 514         cout<<"Opening file "<<filename<<endl;
 515         SetGaliceFile(filename);
 516     }
 517
 518
 519     for (Int_t ievent=fFirstEvent; ievent<=fLastEvent; ievent++) { // event loop
 520         Bool_t isClearEvent = kTRUE;
 521
 522         for(Int_t ch=0; ch<fgkNchambers; ch++){
 523             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 524                 partNumOfTrig[ch][cath]=0;
 525             }
 526             for(Int_t itrack=0; itrack<2; itrack++){
 527                 triggeredDigits[itrack][ch]=-1;
 528             }
 529         }
 530
 531         fRunLoader->GetEvent(ievent);
 532         if (ievent%evtBeforePrint==0) printf("\t Event = %d\n",ievent);
 533
 534         fData->ResetDigits();
 535         fData->ResetTrigger();
 536         fData->ResetRecTriggerTracks();
 537         fData->ResetRecTracks();
 538
 539         fData->SetTreeAddress("RL,RT");
 540         fData->GetRecTriggerTracks();
 541         recTrigTracksArray = fData->RecTriggerTracks();
 542         Int_t nRecTrigTracks = (Int_t) recTrigTracksArray->GetEntriesFast();
 543
 544         fData->SetTreeAddress("D,GLT");
 545         fData->GetDigits();
 546
 547         const AliMUONGeometryTransformer* kGeomTransformer = fMUON->GetGeometryTransformer();
 548
 549         for (Int_t iRecTrigTrack=0; iRecTrigTrack<nRecTrigTracks; iRecTrigTrack++) {
 550             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 551                 digitPerTrack[cath]=0;
 552             }
 553             for(Int_t ch=0; ch<fgkNchambers; ch++){
 554                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 555                     match[ch][cath]=kFALSE;
 556                     slatThatTriggered[ch][cath]=-1;
 557                     for(Int_t loc=0; loc<4; loc++){
 558                         boardThatTriggered[ch][cath][loc]=-1;
 559                     }
 560                 }
 561             }
 562
 563             Bool_t doubleCountTrack = kFALSE;
 564
 565             // reading info from tracks
 566             recTrigTrack = (AliMUONTriggerTrack *)recTrigTracksArray->At(iRecTrigTrack);
 567             Float_t x11 = recTrigTrack->GetX11();// x position (info from non-bending plane)
 568             Float_t y11 = recTrigTrack->GetY11();// y position (info from bending plane)
 569             Float_t thetaX = recTrigTrack->GetThetax();
 570             Float_t thetaY = recTrigTrack->GetThetay();
 571
 572             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);
 573
 574             for(Int_t ch=0; ch<fgkNchambers; ch++) {
 575                 zRealMatch[ch] = zMeanChamber[ch];
 576                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 577                     trigScheme[ch][cath] = 0;
 578                 }
 579             }
 580
 581             for(Int_t ch=0; ch<fgkNchambers; ch++) { // chamber loop
 582                 Int_t currCh = chOrder[ch];
 583                 if(fDebugLevel>=2){
 584                     if(fDebugLevel<3)printf("\tEvent = %i, Track = %i\n", ievent, iRecTrigTrack);
 585                     printf("zMeanChamber[%i] = %.2f\tzRealMatch[0] = %.2f\n",currCh,zMeanChamber[currCh],zRealMatch[0]);
 586                 }
 587
 588                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 589                     correctFactor[cath]=1.;
 590                 }
 591                 // calculate corrections to trigger track theta
 592                 if(ch>=1)correctFactor[kNonBending] = zMeanChamber[0]/zRealMatch[0];// corrects x position
 593                 if(ch>=2)correctFactor[kBending] = (zMeanChamber[2] - zMeanChamber[0]) / (zRealMatch[2] - zRealMatch[0]);// corrects y position
 594
 595                 // searching track intersection with chambers (first approximation)
 596                 Float_t deltaZ = zMeanChamber[currCh] - zMeanChamber[0];
 597                 trackIntersectCh[currCh][0] = zMeanChamber[currCh] * TMath::Tan(thetaX) * correctFactor[kNonBending];// x position (info from non-bending plane)
 598                 trackIntersectCh[currCh][1] = y11 + deltaZ * TMath::Tan(thetaY) * correctFactor[kBending];// y position (info from bending plane)
 599                 Int_t detElemIdFromTrack = DetElemIdFromPos(trackIntersectCh[currCh][0], trackIntersectCh[currCh][1], 11+currCh, 0);
 600                 if(detElemIdFromTrack<0) {
 601                     if(fDebugLevel>1) printf("Warning: trigger track outside trigger chamber\n");
 602                     continue;
 603                 }
 604
 605                 triggeredDigits[1][currCh] = MatchingPad(detElemIdFromTrack, trackIntersectCh[currCh], kGeomTransformer, matchPad, nboard, zRealMatch, y11);
 606
 607                 // deciding if digit matches track
 608                 Bool_t isDiffLocBoard = kFALSE;
 609                 if(fReproduceTrigResponse && ch>2){
 610                     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 611                         if(boardThatTriggered[currCh][cath][0]>=0){
 612                             if(boardThatTriggered[currCh][cath][0]!=boardThatTriggered[currCh-1][cath][0]) isDiffLocBoard = kTRUE;
 613                         }
 614                     }
 615                 }
 616
 617                 if(isDiffLocBoard && fDebugLevel>=1)printf("\tDifferent local board\n");
 618
 619                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 620                     match[currCh][cath] = (matchPad[cath] && !isDiffLocBoard);
 621                     if(!match[currCh][cath]) continue;
 622                     digitPerTrack[cath]++;
 623                     trigScheme[currCh][cath]++;
 624                     slatThatTriggered[currCh][cath] = detElemIdFromTrack;
 625                     for(Int_t loc=0; loc<4; loc++){
 626                         boardThatTriggered[currCh][cath][loc] = nboard[cath][loc];
 627                     }
 628                 }
 629             } // end chamber loop
 630
 631             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 632                 if(digitPerTrack[cath]<3)isClearEvent = kFALSE;
 633                 if(fDebugLevel>=1 && !isClearEvent)printf("Warning: found %i digits for trigger track cathode %i.\nRejecting event\n", digitPerTrack[cath],cath);
 634             }
 635
 636             if(!isClearEvent && !fReproduceTrigResponse) continue;
 637
 638             Int_t commonDigits = 0;
 639             for(Int_t ch=0; ch<fgkNchambers; ch++){
 640                 if(triggeredDigits[1][ch]==triggeredDigits[0][ch]) commonDigits++; // Compare with previous track
 641                 triggeredDigits[0][ch] = triggeredDigits[1][ch]; // Store this track parameters for comparison with next one
 642             }
 643             if(commonDigits>=2){
 644                 doubleCountTrack=kTRUE;
 645             }
 646
 647             if(!doubleCountTrack || fReproduceTrigResponse){
 648                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 649                     Int_t is44 = 1;
 650                     Bool_t goodForSlatEff = kTRUE;
 651                     Bool_t goodForBoardEff = kTRUE;
 652                     Int_t ineffSlat = -1;
 653                     Int_t ineffDetElId = -1;
 654                     Int_t firstSlat = slatThatTriggered[0][cath]%100;
 655                     if(firstSlat<0) firstSlat=slatThatTriggered[1][cath]%100;
 656                     Int_t firstBoard = boardThatTriggered[0][kBending][0];
 657                     if(firstBoard<0) firstBoard=boardThatTriggered[1][kBending][0];
 658                     for(Int_t ch=0; ch<fgkNchambers; ch++){
 659                         Bool_t isCurrChIneff = kFALSE;
 660                         is44 *= trigScheme[ch][cath];
 661                         Int_t currSlat = slatThatTriggered[ch][cath]%100;
 662                         if(currSlat<0){
 663                             ineffDetElId = DetElemIdFromPos(trackIntersectCh[ch][0], trackIntersectCh[ch][1], 11+ch, cath);
 664                             currSlat = ineffDetElId%100;
 665                             ineffSlat = currSlat;
 666                             isCurrChIneff = kTRUE;
 667                         }
 668                         if(currSlat!=firstSlat)goodForSlatEff=kFALSE;
 669                         Bool_t atLeastOneLoc=kFALSE;
 670                         if(isCurrChIneff) LocalBoardFromPos(trackIntersectCh[ch][0], trackIntersectCh[ch][1], ineffDetElId, cath, ineffBoard);
 671                         for(Int_t loc=0; loc<4; loc++){
 672                             Int_t currBoard = boardThatTriggered[ch][cath][loc];
 673                             if(isCurrChIneff) currBoard = ineffBoard[loc];
 674                             if(currBoard==firstBoard){
 675                                 atLeastOneLoc=kTRUE;
 676                                 break;
 677                             }
 678                         }
 679                         if(!atLeastOneLoc)goodForBoardEff=kFALSE;
 680                     } // end chamber loop
 681
 682                     for(Int_t ch=0; ch<fgkNchambers; ch++){
 683                         if(match[ch][cath])partNumOfTrig[ch][cath]++;
 684                     }
 685
 686                     // Trigger 4/4
 687                     if(is44==1){
 688                         fTrigger44[cath]++;
 689                         if(fDebugLevel>=1)printf("Trigger44[%i] = %i\n",cath,fTrigger44[cath]);
 690                         if(goodForSlatEff){
 691                             for(Int_t ch=0; ch<fgkNchambers; ch++){
 692                                 fHitPerSlat[ch][cath][firstSlat]++;
 693                                 if(fDebugLevel>=1)printf("Slat that triggered = %i\n",slatThatTriggered[ch][cath]);
 694                                 if(goodForBoardEff && firstBoard>0){
 695                                     fHitPerBoard[ch][cath][firstBoard-1]++;
 696                                     if(fDebugLevel>=1)printf("Board that triggered = %i\n",firstBoard);
 697                                 }
 698                                 else if(fDebugLevel>=1)printf("Event = %i, Track = %i: Particle crossed different boards: rejected!\n",ievent,iRecTrigTrack);
 699                             }
 700                         }
 701                         else printf("Event = %i, Track = %i: Particle crossed different slats: rejected!\n",ievent,iRecTrigTrack);
 702                     }
 703
 704                     // Trigger 3/4
 705                     if(ineffDetElId>0){
 706                         Int_t ineffCh = ineffDetElId/100-11;
 707                         fTrigger34[ineffCh][cath]++;
 708                         if(fDebugLevel>=1) printf("Trigger34[%i][%i] = %i\n",ineffCh,cath,fTrigger34[ineffCh][cath]);
 709                         if(goodForSlatEff){
 710                             if(fDebugLevel>=1) printf("Slat non efficient = %i\n",ineffDetElId);
 711                             fInefficientSlat[ineffCh][cath][ineffSlat]++;
 712
 713                             if(goodForBoardEff && firstBoard>0){
 714                                 if(fDebugLevel>=1) printf("Board non efficient = %i\n",firstBoard);
 715                                 fInefficientBoard[ineffCh][cath][firstBoard-1]++;
 716                             }
 717                             else if(fDebugLevel>=1) printf("Event = %i, Track = %i: Particle crossed different boards: rejected!\n",ievent,iRecTrigTrack);
 718                         }
 719                         else printf("Event %i, Track = %i: Particle crossed different slats: rejected!\n",ievent,iRecTrigTrack);
 720                     }
 721                 } // end loop on cathodes
 722             }
 723             else if(doubleCountTrack){
 724                 if(fDebugLevel<=1)printf("\n\tEvent = %i, Track = %i: ", ievent,iRecTrigTrack);
 725                 printf("Double Count Track: %i similar to %i. Track rejected!\n",iRecTrigTrack, iRecTrigTrack-1);
 726             }
 727         } // end trigger tracks loop
 728         if(nRecTrigTracks<=0) continue;
 729
 730         for(Int_t ch=0; ch<fgkNchambers; ch++){
 731             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 732                 totNumOfTrig[ch][cath] += partNumOfTrig[ch][cath];
 733                 if(partNumOfTrig[ch][cath]>0)atLeast1MuPerEv[ch][cath]++;
 734             }
 735         }
 736
 737         if(fPrintInfo){
 738             //Global trigger
 739             globalTrigger = fData->GlobalTrigger();
 740             Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
 741
 742             for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
 743                 gloTrg = (AliMUONGlobalTrigger*)globalTrigger->At(iglobal);
 744             }
 745             PrintTrigger(gloTrg);
 746             InfoDigit();
 747             cout<<"\n"<<endl;
 748         }
 749     }// end event loop
 750     if(fFirstRun>=0)CleanGalice();
 751     } //end loop over run
 752
 753     // Write output data
 754     WriteEfficiencyMap(outputDir);
 755
 756     WriteOutput(outputDir, totNumOfTrig, atLeast1MuPerEv);
 757 }
 758
 759 //_____________________________________________________________________________
 760 void AliMUONTriggerChamberEff::WriteOutput(const char* outputDir, Int_t totNumOfTrig[4][2], Int_t atLeast1MuPerEv[4][2])
 761 {
 762     //
 763     /// Writes information on calculated efficiency.
 764     /// It writes: triggerChamberEff.root file containing efficiency histograms.
 765     ///
 766     /// In addition a text file triggerChamberEff.out is created,
 767     /// with further informations on efficiencies.
 768     //
 769
 770     char *cathodeName[fgkNcathodes]={"Bending plane", "Non-Bending plane"};
 771     char *cathCode[fgkNcathodes] = {"bendPlane", "nonBendPlane"};
 772
 773     char outFileName[100];
 774     sprintf(outFileName, "%s/triggerChamberEff.out",outputDir);
 775     FILE *outfile = fopen(outFileName, "w");
 776     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 777         fprintf(outfile,"%s:\n",cathodeName[cath]);
 778         for(Int_t ch=0; ch<fgkNchambers; ch++){
 779             fprintf(outfile,"Total number of muon triggers chamber 1%i = %i\n",ch+1,totNumOfTrig[ch][cath]);
 780         }
 781         fprintf(outfile,"\n");
 782     }
 783     fprintf(outfile,"\n");
 784     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 785         fprintf(outfile,"%s:\n",cathodeName[cath]);
 786         for(Int_t ch=0; ch<fgkNchambers; ch++){
 787             fprintf(outfile,"At least 1 muon triggered chamber 1%i = %i\n",ch+1, atLeast1MuPerEv[ch][cath]);
 788         }
 789         fprintf(outfile,"\n");
 790     }
 791     fprintf(outfile,"\n\n");
 792     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 793         fprintf(outfile,"%s:\n",cathodeName[cath]);
 794         fprintf(outfile,"Number of triggers where all chambers counted = %i\n",fTrigger44[cath]);
 795         fprintf(outfile,"\n");
 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             fprintf(outfile,"Number of triggers where chamber 1%i did not count = %i\n",ch+1,fTrigger34[ch][cath]);
 802         }
 803         fprintf(outfile,"\n");
 804     }
 805
 806     fprintf(outfile,"\n");
 807     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 808         fprintf(outfile,"%s:\n",cathodeName[cath]);
 809         for(Int_t ch=0; ch<fgkNchambers; ch++){
 810             Int_t sumIneff = 0, sumHits = 0;
 811             fprintf(outfile,"\n Chamber %1i\n", ch+1);
 812             for(Int_t slat=0; slat<fgkNslats; slat++){
 813                 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]);
 814                 sumIneff += fInefficientSlat[ch][cath][slat];
 815                 sumHits += fHitPerSlat[ch][cath][slat];
 816             }
 817             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);
 818         }
 819         fprintf(outfile,"\n");
 820     }
 821     fclose(outfile);
 822
 823     sprintf(outFileName, "%s/triggerChamberEff.root",outputDir);
 824     TFile *outputHistoFile = new TFile(outFileName,"RECREATE");
 825     TDirectory *dir = gDirectory;
 826
 827     enum {kSlatIn11, kSlatIn12, kSlatIn13, kSlatIn14, kChamberEff};
 828     char *yAxisTitle = "trigger efficiency (a.u.)";
 829     char *xAxisTitle = "chamber";
 830
 831     TH1F *histo[fgkNcathodes][fgkNchambers+1];
 832     TH1F *histoBoard[fgkNcathodes][fgkNchambers];
 833
 834     char histoName[30];
 835     char histoTitle[90];
 836
 837     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 838         for(Int_t ch=0; ch<fgkNchambers+1; ch++){
 839             if(ch==kChamberEff){
 840                 sprintf(histoName, "%sChamberEff", cathCode[cath]);
 841                 sprintf(histoTitle, "Chamber efficiency %s", cathCode[cath]);
 842                 histo[cath][ch] = new TH1F(histoName, histoTitle, fgkNchambers, 11-0.5, 15-0.5);
 843                 histo[cath][ch]->SetXTitle(xAxisTitle);
 844                 histo[cath][ch]->SetYTitle(yAxisTitle);
 845                 histo[cath][ch]->GetXaxis()->SetNdivisions(fgkNchambers);
 846             }
 847             else {
 848                 sprintf(histoName, "%sSlatEffChamber%i", cathCode[cath], 11+ch);
 849                 sprintf(histoTitle, "Chamber %i: slat efficiency %s", 11+ch, cathCode[cath]);
 850                 histo[cath][ch] = new TH1F(histoName, histoTitle, fgkNslats, 0-0.5, fgkNslats-0.5);
 851                 histo[cath][ch]->SetXTitle("slat");
 852                 histo[cath][ch]->SetYTitle(yAxisTitle);
 853                 histo[cath][ch]->GetXaxis()->SetNdivisions(fgkNslats);
 854
 855                 sprintf(histoName, "%sBoardEffChamber%i", cathCode[cath], 11+ch);
 856                 sprintf(histoTitle, "Chamber %i: board efficiency %s", 11+ch, cathCode[cath]);
 857                 histoBoard[cath][ch] = new TH1F(histoName, histoTitle, fgkNboards, 1-0.5, fgkNboards+1.-0.5);
 858                 histoBoard[cath][ch]->SetXTitle("boards");
 859                 histoBoard[cath][ch]->SetYTitle(yAxisTitle);
 860                 histoBoard[cath][ch]->GetXaxis()->SetNdivisions(fgkNboards);
 861             }
 862         }
 863     }
 864
 865     Float_t efficiency, efficiencyError;
 866
 867     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 868         for(Int_t ch=0; ch<fgkNchambers; ch++){
 869             for(Int_t slat=0; slat<fgkNslats; slat++){
 870                 CalculateEfficiency(fHitPerSlat[ch][cath][slat], fHitPerSlat[ch][cath][slat]+fInefficientSlat[ch][cath][slat], efficiency, efficiencyError, kFALSE);
 871                 histo[cath][ch]->SetBinContent(slat+1, efficiency);
 872                 histo[cath][ch]->SetBinError(slat+1, efficiencyError);
 873             }
 874             CalculateEfficiency(fTrigger44[cath], fTrigger34[ch][cath]+fTrigger44[cath], efficiency, efficiencyError, kFALSE);
 875             histo[cath][kChamberEff]->SetBinContent(ch+1, efficiency);
 876             histo[cath][kChamberEff]->SetBinError(ch+1, efficiencyError);
 877
 878             for(Int_t board=0; board<fgkNboards; board++){
 879                 CalculateEfficiency(fHitPerBoard[ch][cath][board], fHitPerBoard[ch][cath][board]+fInefficientBoard[ch][cath][board], efficiency, efficiencyError, kFALSE);
 880                 histoBoard[cath][ch]->SetBinContent(board+1, efficiency);
 881                 histoBoard[cath][ch]->SetBinError(board+1, efficiencyError);
 882             }
 883         }
 884     }
 885
 886 // write all histos
 887     outputHistoFile->cd();
 888     dir->GetList()->Write();
 889     outputHistoFile->Close();
 890 }
 891
 892
 893 //_____________________________________________________________________________
 894 void AliMUONTriggerChamberEff::WriteEfficiencyMap(const char* outputDir)
 895 {
 896     //
 897     /// Writes the calculated efficiency in the text file efficiencyCells.dat
 898     ///
 899     /// The file can be further put in $ALICE_ROOT/MUON/data
 900     /// and used to run simulations with measured trigger chamber efficiencies.
 901     //
 902
 903     Int_t effOutWidth=4;
 904
 905     Float_t efficiency, efficiencyError;
 906
 907     Int_t aCapo[] = {16, 38, 60, 76, 92, 108, 117, 133, 155, 177, 193, 209, 225, 234};
 908
 909     char filename[70];
 910     sprintf(filename, "%s/efficiencyCells.dat", outputDir);
 911     ofstream outFile(filename);
 912     outFile << "localBoards" << endl;
 913     for(Int_t ch=0; ch<fgkNchambers; ch++){
 914         //Print information
 915         outFile << "\n\ndetElemId:\t" << 11+ch;
 916         outFile << "00";
 917         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 918             outFile << "\n cathode:\t" << cath << endl;
 919             Int_t currLine=0;
 920             for(Int_t board=0; board<fgkNboards; board++){
 921
 922                 if(board==aCapo[currLine]){
 923                     outFile << endl;
 924                     currLine++;
 925                 }
 926                 CalculateEfficiency(fHitPerBoard[ch][cath][board], fHitPerBoard[ch][cath][board]+fInefficientBoard[ch][cath][board], efficiency, efficiencyError, kFALSE);
 927                 outFile << " " << setw(effOutWidth) << efficiency;
 928             }// loop on boards
 929             outFile << endl;
 930         }// loop on cathodes
 931     }// loop on chambers
 932 }
 933