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