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 "AliMUONVDigit.h"
  27 #include "AliMUONConstants.h"
  28 #include "AliMUONTriggerTrack.h"
  29 #include "AliMUONDigitMaker.h"
  30 #include "AliMUONLocalTrigger.h"
  31 #include "AliMUONGeometryTransformer.h"
  32
  33 #include "AliMUONTrack.h"
  34 #include "AliMUONTrackParam.h"
  35 #include "AliMUONTrackExtrap.h"
  36
  37 #include "AliMUONDigitStoreV1.h"
  38 #include "AliMUONVDigitStore.h"
  39 #include "AliMUONVTriggerStore.h"
  40 #include "AliMUONVTriggerTrackStore.h"
  41 #include "AliMUONVTrackStore.h"
  42
  43 #include "AliMpVSegmentation.h"
  44 #include "AliMpSegmentation.h"
  45 #include "AliMpPad.h"
  46 #include "AliMpDEIterator.h"
  47 #include "AliMpPlaneType.h"
  48 #include "AliMpDEManager.h"
  49
  50 #include "AliLog.h"
  51
  52 #include <Riostream.h>
  53 #include <TFile.h>
  54 #include <TH1F.h>
  55 #include <TMath.h>
  56
  57 #include <TSeqCollection.h>
  58 #include <TTree.h>
  59 #include <TROOT.h>
  60
  61
  62 /// \cond CLASSIMP
  63 ClassImp(AliMUONTriggerChamberEff)
  64 /// \endcond
  65
  66
  67 //_____________________________________________________________________________
  68 AliMUONTriggerChamberEff::AliMUONTriggerChamberEff()
  69 : TObject(),
  70   fTransformer(0x0),
  71   fDigitMaker(0x0),
  72   fReproduceTrigResponse(kFALSE),
  73   fPrintInfo(kFALSE),
  74   fWriteOnESD(kFALSE),
  75   fDebugLevel(0),
  76   fkMaxDistance(99999.)
  77 {
  78 /// Standard constructor
  79     ResetArrays();
  80 }
  81
  82
  83 //_____________________________________________________________________________
  84 AliMUONTriggerChamberEff::AliMUONTriggerChamberEff(const AliMUONGeometryTransformer* transformer,
  85                                                    const AliMUONDigitMaker* digitMaker,
  86                                                    Bool_t writeOnESD)
  87 : TObject(),
  88   fTransformer(transformer),
  89   fDigitMaker(digitMaker),
  90   fReproduceTrigResponse(kFALSE),
  91   fPrintInfo(kFALSE),
  92   fWriteOnESD(writeOnESD),
  93   fDebugLevel(0),
  94   fkMaxDistance(99999.)
  95 {
  96 /// Standard constructor
  97     ResetArrays();
  98 }
  99
 100
 101 //_____________________________________________________________________________
 102 AliMUONTriggerChamberEff::~AliMUONTriggerChamberEff()
 103 {
 104 /// Destructor
 105     Bool_t writeOnESD=fWriteOnESD;
 106     fWriteOnESD=kFALSE;
 107     if(writeOnESD) SaveInESDFile();
 108 }
 109
 110
 111 //_____________________________________________________________________________
 112 AliMUONTriggerChamberEff::AliMUONTriggerChamberEff(const AliMUONTriggerChamberEff& other)
 113     :TObject(other),
 114      fTransformer(0x0),
 115      fDigitMaker(0x0),
 116      fReproduceTrigResponse(other.fReproduceTrigResponse),
 117      fPrintInfo(other.fPrintInfo),
 118      fWriteOnESD(other.fWriteOnESD),
 119      fDebugLevel(other.fDebugLevel),
 120      fkMaxDistance(other.fkMaxDistance)
 121 {
 122 /// copy constructor
 123
 124     for(Int_t ch=0; ch<fgkNchambers; ch++){
 125         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 126             fTrigger34[ch][cath] = other.fTrigger34[ch][cath];
 127             fTrigger44[cath] = other.fTrigger44[cath];
 128             for(Int_t slat=0; slat<fgkNslats; slat++){
 129                 fInefficientSlat[ch][cath][slat] = other.fInefficientSlat[ch][cath][slat];
 130                 fHitPerSlat[ch][cath][slat] = other.fHitPerSlat[ch][cath][slat];
 131             }
 132             for(Int_t board=0; board<fgkNboards; board++){
 133                 fInefficientBoard[ch][cath][board] = other.fInefficientBoard[ch][cath][board];
 134                 fHitPerBoard[ch][cath][board] = other.fHitPerBoard[ch][cath][board];
 135             }
 136         }
 137     }
 138 }
 139
 140
 141 //_____________________________________________________________________________
 142 AliMUONTriggerChamberEff& AliMUONTriggerChamberEff::operator=(const AliMUONTriggerChamberEff& other)
 143 {
 144     /// Asignment operator
 145     // check assignement to self
 146     if (this == &other)
 147         return *this;
 148
 149     // base class assignement
 150     TObject::operator=(other);
 151
 152     fTransformer = 0x0;
 153     fDigitMaker = 0x0;
 154     fReproduceTrigResponse = other.fReproduceTrigResponse;
 155     fPrintInfo = other.fPrintInfo;
 156     fWriteOnESD = other.fWriteOnESD;
 157     fDebugLevel = other.fDebugLevel;
 158     //fkMaxDistance = other.fkMaxDistance;
 159     return *this;
 160 }
 161
 162
 163 //_____________________________________________________________________________
 164 void AliMUONTriggerChamberEff::ResetArrays()
 165 {
 166     //
 167     /// Sets the data member counters to 0.
 168     //
 169
 170     for(Int_t ch=0; ch<fgkNchambers; ch++){
 171         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 172             fTrigger34[ch][cath] = 0;
 173             fTrigger44[cath] = 0;
 174             for(Int_t slat=0; slat<fgkNslats; slat++){
 175                 fInefficientSlat[ch][cath][slat] = 0;
 176                 fHitPerSlat[ch][cath][slat] = 0;
 177             }
 178             for(Int_t board=0; board<fgkNboards; board++){
 179                 fInefficientBoard[ch][cath][board] = 0;
 180                 fHitPerBoard[ch][cath][board] = 0;
 181             }
 182         }
 183     }
 184 }
 185
 186
 187 //______________________________________________________________________________
 188 Bool_t
 189 AliMUONTriggerChamberEff::TriggerDigits(const AliMUONVTriggerStore& triggerStore,
 190                                         AliMUONVDigitStore& digitStore) const
 191 {
 192     //
 193     /// make (S)Digit for trigger
 194     //
 195     digitStore.Clear();
 196
 197     AliMUONLocalTrigger* locTrg;
 198     TIter next(triggerStore.CreateLocalIterator());
 199
 200     while ( ( locTrg = static_cast<AliMUONLocalTrigger*>(next()) ) )
 201     {
 202         if (locTrg->IsNull()) continue;
 203
 204         TArrayS xyPattern[2];
 205         locTrg->GetXPattern(xyPattern[0]);
 206         locTrg->GetYPattern(xyPattern[1]);
 207
 208         Int_t nBoard = locTrg->LoCircuit();
 209         fDigitMaker->TriggerDigits(nBoard, xyPattern, digitStore);
 210     }
 211     return kTRUE;
 212 }
 213
 214
 215 //_____________________________________________________________________________
 216 void AliMUONTriggerChamberEff::InfoDigit(AliMUONVDigitStore& digitStore)
 217 {
 218     //
 219     /// Prints information on digits (for debugging)
 220     //
 221     TIter next(digitStore.CreateIterator());
 222     AliMUONVDigit* mDigit=0x0;
 223
 224     while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
 225     {
 226         mDigit->Print();
 227     } // end digit loop
 228     printf("\n");
 229 }
 230
 231
 232 //_____________________________________________________________________________
 233 Int_t AliMUONTriggerChamberEff::MatchingPad(AliMUONVDigitStore& digitStore, Int_t &detElemId,
 234                                             Float_t coor[2], Bool_t isMatch[fgkNcathodes],
 235                                             Int_t nboard[fgkNcathodes][4],
 236                                             Float_t zRealMatch[fgkNchambers], Float_t y11)
 237 {
 238     //
 239     /// Check slat and board number of digit matching track
 240     //
 241
 242     enum {kBending, kNonBending};
 243
 244     Float_t minMatchDist[fgkNcathodes];
 245     Int_t padsInCheckArea[fgkNcathodes];
 246
 247     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 248         isMatch[cath] = kFALSE;
 249         minMatchDist[cath] = fkMaxDistance/10.;
 250         padsInCheckArea[cath] = 0;
 251     }
 252     Int_t iChamber = AliMpDEManager::GetChamberId(detElemId);
 253     Int_t ch = iChamber-10;
 254     Float_t oldDeltaZ = AliMUONConstants::DefaultChamberZ(iChamber) - AliMUONConstants::DefaultChamberZ(10);
 255     Float_t y = coor[1];
 256     Int_t iSlat = detElemId%100;
 257     Int_t trigDigitBendPlane = -1;
 258     Int_t foundDetElemId = detElemId;
 259     Float_t foundZmatch=999.;
 260     Float_t yCoorAtPadZ=999.;
 261
 262     TIter next(digitStore.CreateIterator());
 263     AliMUONVDigit* mDigit;
 264     Int_t idigit=0;
 265
 266     while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
 267     {
 268         idigit++;
 269         Int_t currDetElemId = mDigit->DetElemId();
 270         Int_t currCh = AliMpDEManager::GetChamberId(currDetElemId);
 271         if(currCh!=iChamber) continue;
 272         Int_t currSlat = currDetElemId%100;
 273         Int_t slatDiff = TMath::Abs(currSlat-iSlat);
 274         if(slatDiff>1 && slatDiff<17) continue; // Check neighbour slats
 275         Int_t cathode = mDigit->Cathode();
 276         Int_t ix = mDigit->PadX();
 277         Int_t iy = mDigit->PadY();
 278         Float_t xpad, ypad, zpad;
 279         const AliMpVSegmentation* seg = AliMpSegmentation::Instance()
 280             ->GetMpSegmentation(currDetElemId,AliMp::GetCathodType(cathode));
 281
 282         AliMpPad pad = seg->PadByIndices(AliMpIntPair(ix,iy),kTRUE);
 283         Float_t xlocal1 = pad.Position().X();
 284         Float_t ylocal1 = pad.Position().Y();
 285         Float_t dpx = pad.Dimensions().X();
 286         Float_t dpy = pad.Dimensions().Y();
 287         fTransformer->Local2Global(currDetElemId, xlocal1, ylocal1, 0, xpad, ypad, zpad);
 288         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]);
 289         // searching track intersection with chambers (second approximation)
 290         if(ch%2==1){
 291             //if(iChamber%2==1){
 292             Float_t deltaZ = zpad - zRealMatch[0];
 293             y = (coor[1]-y11)*deltaZ/oldDeltaZ + y11;
 294             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);
 295         }
 296         Float_t matchDist = PadMatchTrack(xpad, ypad, dpx, dpy, coor[0], y, ch);
 297         if(matchDist<fkMaxDistance/2.) padsInCheckArea[cathode]++;
 298         if(matchDist>minMatchDist[cathode])continue;
 299         isMatch[cathode] = kTRUE;
 300         minMatchDist[cathode] = matchDist;
 301         foundDetElemId = currDetElemId;
 302         foundZmatch=zpad;
 303         yCoorAtPadZ=y;
 304         if(cathode==kBending) trigDigitBendPlane = idigit;
 305         for (Int_t loc=0; loc<pad.GetNofLocations(); loc++){
 306             AliMpIntPair location = pad.GetLocation(loc);
 307             nboard[cathode][loc] = location.GetFirst();
 308         }
 309         for(Int_t loc=pad.GetNofLocations(); loc<4; loc++){
 310             nboard[cathode][loc]=-1;
 311         }
 312     }
 313
 314     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 315         if(padsInCheckArea[cath]>2) {
 316             if(fDebugLevel>=1) printf("padsInCheckArea[%i] = %i\n",cath,padsInCheckArea[cath]);
 317             return -500;
 318         }
 319     }
 320
 321     if(isMatch[kBending] || isMatch[kNonBending]){
 322         detElemId = foundDetElemId;
 323         zRealMatch[ch] = foundZmatch;
 324         coor[1] = yCoorAtPadZ;
 325         if(fDebugLevel>2){
 326             Int_t whichCathode=kBending;
 327             if(!isMatch[kBending])whichCathode=kNonBending;
 328         }
 329     }
 330     return trigDigitBendPlane;
 331 }
 332
 333 //_____________________________________________________________________________
 334 Float_t AliMUONTriggerChamberEff::PadMatchTrack(Float_t xPad, Float_t yPad,
 335                                                 Float_t dpx, Float_t dpy,
 336                                                 Float_t xTrackAtPad, Float_t yTrackAtPad,
 337                                                 Int_t chamber)
 338 {
 339     //
 340     /// Decides if the digit belongs to the trigger track.
 341     //
 342
 343     Float_t maxDist = 2.;//3. // cm
 344     Float_t maxDistCheckArea = 6.; // cm
 345
 346     Float_t matchDist = fkMaxDistance;
 347
 348     Float_t deltaX = TMath::Abs(xPad-xTrackAtPad)-dpx;
 349     Float_t deltaY = TMath::Abs(yPad-yTrackAtPad)-dpy;
 350     Float_t maxDistX = maxDist;
 351     Float_t maxDistY = maxDist;
 352
 353     if(fReproduceTrigResponse){
 354         maxDistX = dpx;
 355         maxDistY = dpy;
 356         deltaX = TMath::Abs(xPad-xTrackAtPad);
 357         deltaY = TMath::Abs(yPad-yTrackAtPad);
 358         if(dpx<dpy && chamber>=2) maxDistX = 3.*dpx;// Non-bending plane: check the +- 1 strip between stations
 359         if(dpy<dpx && chamber%2) maxDistY = 3.*dpy;// bending plane: check the +- 1 strip between planes in the same station
 360     }
 361
 362     if(deltaX<=maxDistX && deltaY<=maxDistY) matchDist = TMath::Max(deltaX, deltaY);
 363     else if(deltaX<=maxDistCheckArea && deltaY<=maxDistCheckArea) matchDist = fkMaxDistance/5.;
 364     return matchDist;
 365 }
 366
 367
 368 //_____________________________________________________________________________
 369 void AliMUONTriggerChamberEff::CalculateEfficiency(Int_t trigger44, Int_t trigger34,
 370                                                    Float_t &efficiency, Float_t &error,
 371                                                    Bool_t failuresAsInput)
 372 {
 373     //
 374     /// Returns the efficiency.
 375     //
 376
 377     efficiency=-9.;
 378     error=0.;
 379     if(trigger34>0){
 380         efficiency=(Double_t)trigger44/((Double_t)trigger34);
 381         if(failuresAsInput)efficiency=1.-(Double_t)trigger44/((Double_t)trigger34);
 382     }
 383     Double_t q = TMath::Abs(1-efficiency);
 384     if(efficiency<0)error=0.0;
 385     else error = TMath::Sqrt(efficiency*q/((Double_t)trigger34));
 386 }
 387
 388
 389 //_____________________________________________________________________________
 390 Int_t AliMUONTriggerChamberEff::DetElemIdFromPos(Float_t x, Float_t y,
 391                                                  Int_t chamber, Int_t cathode)
 392 {
 393     //
 394     /// Given the (x,y) position in the chamber,
 395     /// it returns the corresponding slat
 396     //
 397
 398     Int_t resultingDetElemId = -1;
 399     AliMpDEIterator it;
 400     Float_t minDist = 999.;
 401     for ( it.First(chamber-1); ! it.IsDone(); it.Next() ){
 402         Int_t detElemId = it.CurrentDEId();
 403         Int_t ich = detElemId/100-10;
 404         Float_t tolerance=0.2*((Float_t)ich);
 405         Float_t currDist=9999.;
 406
 407         const AliMpVSegmentation* seg =
 408             AliMpSegmentation::Instance()
 409             ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 410         if (!seg) continue;
 411
 412         Float_t deltax = seg->Dimensions().X();
 413         Float_t deltay = seg->Dimensions().Y();
 414         Float_t xlocal1 =  -deltax;
 415         Float_t ylocal1 =  -deltay;
 416         Float_t xlocal2 =  +deltax;
 417         Float_t ylocal2 =  +deltay;
 418         Float_t xg01, yg01, zg1, xg02, yg02, zg2;
 419         fTransformer->Local2Global(detElemId, xlocal1, ylocal1, 0, xg01, yg01, zg1);
 420         fTransformer->Local2Global(detElemId, xlocal2, ylocal2, 0, xg02, yg02, zg2);
 421
 422         Float_t xg1 = xg01, xg2 = xg02, yg1 = yg01, yg2 = yg02;
 423
 424         if(xg01>xg02){
 425             xg1 = xg02;
 426             xg2 = xg01;
 427         }
 428         if(yg01>yg02){
 429             yg1 = yg02;
 430             yg2 = yg01;
 431         }
 432
 433         if(x>=xg1-tolerance && x<=xg2+tolerance && y>=yg1-tolerance && y<=yg2+tolerance){ // takes into account errors in extrapolation
 434             if(y<yg1) currDist = yg1-y;
 435             else if(y>yg2) currDist = y-yg2;
 436             if(currDist<minDist) {
 437                 resultingDetElemId = detElemId;
 438                 minDist=currDist;
 439                 continue;
 440             }
 441             resultingDetElemId = detElemId;
 442             break;
 443         }
 444     } // loop on detElemId
 445     return resultingDetElemId;
 446 }
 447
 448
 449 //_____________________________________________________________________________
 450 void AliMUONTriggerChamberEff::LocalBoardFromPos(Float_t x, Float_t y,
 451                                                  Int_t detElemId, Int_t cathode,
 452                                                  Int_t localBoard[4])
 453 {
 454     //
 455     /// Given the (x,y) position in the chamber,
 456     /// it returns the corresponding local board
 457     //
 458
 459     for(Int_t loc=0; loc<4; loc++){
 460         localBoard[loc]=-1;
 461     }
 462     Float_t xl, yl, zl;
 463     fTransformer->Global2Local(detElemId, x, y, 0, xl, yl, zl);
 464     TVector2 pos(xl,yl);
 465     const AliMpVSegmentation* seg =
 466         AliMpSegmentation::Instance()
 467           ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 468     if (seg){
 469         AliMpPad pad = seg->PadByPosition(pos,kFALSE);
 470         for (Int_t loc=0; loc<pad.GetNofLocations(); loc++){
 471             AliMpIntPair location = pad.GetLocation(loc);
 472             localBoard[loc] = location.GetFirst();
 473         }
 474     }
 475 }
 476
 477
 478 //_____________________________________________________________________________
 479 void AliMUONTriggerChamberEff::EventChamberEff(const AliMUONVTriggerStore& triggerStore,
 480                                                const AliMUONVTriggerTrackStore& trigTrackStore,
 481                                                const AliMUONVTrackStore& trackStore)
 482 {
 483     //
 484     /// Main method.
 485     /// It loops over the the trigger rec. tracks in the event.
 486     /// Then it search for matching digits around the track.
 487     /// Finally it calculates the efficiency for each trigger board.
 488     /// Files with calculated efficiency are placed in the user defined outputDir.
 489     //
 490
 491     if(!fTransformer || ! fDigitMaker) {
 492         AliError(Form("AliMUONGeometryTransformer or AliMUONDigitMaker not properly initialized!!"));
 493         return;
 494     }
 495
 496     enum {kBending, kNonBending};
 497     Float_t rad2deg = 180./TMath::Pi();
 498
 499     Int_t chOrder[fgkNchambers] = {0,2,1,3};
 500     Float_t zRealMatch[fgkNchambers] = {0.0};
 501     Float_t correctFactor[fgkNcathodes] = {1.};
 502
 503     Bool_t match[fgkNchambers][fgkNcathodes] = {{kFALSE}};
 504     Bool_t matchPad[fgkNcathodes]={kFALSE};
 505
 506
 507     Float_t zMeanChamber[fgkNchambers];
 508     for(Int_t ch=0; ch<fgkNchambers; ch++){
 509         zMeanChamber[ch] = AliMUONConstants::DefaultChamberZ(10+ch);
 510     }
 511
 512     Int_t digitPerTrack[fgkNcathodes] = {0};
 513
 514     Float_t trackIntersectCh[fgkNchambers][2]={{0.0}};
 515
 516     Int_t triggeredDigits[2][fgkNchambers] = {{-1}};
 517
 518     Int_t trigScheme[fgkNchambers][fgkNcathodes]={{0}};
 519     Int_t slatThatTriggered[fgkNchambers][fgkNcathodes]={{-1}};
 520     Int_t boardThatTriggered[fgkNchambers][fgkNcathodes][4]={{{-1}}};
 521     Int_t nboard[fgkNcathodes][4]={{-1}};
 522     Int_t ineffBoard[4]={-1};
 523
 524     AliMUONDigitStoreV1 digitStore;
 525     TriggerDigits(triggerStore,digitStore);
 526
 527     for(Int_t ch=0; ch<fgkNchambers; ch++){
 528         for(Int_t itrack=0; itrack<2; itrack++){
 529             triggeredDigits[itrack][ch]=-1;
 530         }
 531     }
 532
 533     AliMUONTriggerTrack *recTrigTrack = 0x0;
 534
 535     TIter next(trigTrackStore.CreateIterator());
 536
 537     while ( ( recTrigTrack = static_cast<AliMUONTriggerTrack*>(next()) ) )
 538     {
 539         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 540             digitPerTrack[cath]=0;
 541         }
 542         for(Int_t ch=0; ch<fgkNchambers; ch++){
 543             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 544                 match[ch][cath]=kFALSE;
 545                 slatThatTriggered[ch][cath]=-1;
 546                 for(Int_t loc=0; loc<4; loc++){
 547                     boardThatTriggered[ch][cath][loc]=-1;
 548                 }
 549             }
 550         }
 551
 552         Bool_t isClearEvent = kTRUE;
 553         Bool_t doubleCountTrack = kFALSE;
 554
 555         if(!IsCleanTrack(recTrigTrack, trackStore)) {
 556             if(fDebugLevel>=1) printf("\tTrack %p (%f, %f) don't match tracker track: rejected!\n",recTrigTrack,recTrigTrack->GetX11(),recTrigTrack->GetY11());
 557             continue;
 558         }
 559
 560         Float_t x11 = recTrigTrack->GetX11();// x position (info from non-bending plane)
 561         Float_t y11 = recTrigTrack->GetY11();// y position (info from bending plane)
 562         Float_t thetaX = recTrigTrack->GetThetax();
 563         Float_t thetaY = recTrigTrack->GetThetay();
 564
 565         if(fDebugLevel>=3) printf("\tTrack = %p\npos from track: (x,y) = (%f, %f), (thetaX, thetaY) = (%f, %f)\n",recTrigTrack,x11,y11,thetaX*rad2deg,thetaY*rad2deg);
 566
 567         for(Int_t ch=0; ch<fgkNchambers; ch++) {
 568             zRealMatch[ch] = zMeanChamber[ch];
 569             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 570                 trigScheme[ch][cath] = 0;
 571             }
 572         }
 573
 574         for(Int_t ch=0; ch<fgkNchambers; ch++) { // chamber loop
 575             Int_t currCh = chOrder[ch];
 576             if(fDebugLevel>=2)
 577                 printf("zMeanChamber[%i] = %.2f\tzRealMatch[0] = %.2f\n",currCh,zMeanChamber[currCh],zRealMatch[0]);
 578
 579             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 580                 correctFactor[cath]=1.;
 581             }
 582             // calculate corrections to trigger track theta
 583             if(ch>=1)correctFactor[kNonBending] = zMeanChamber[0]/zRealMatch[0];// corrects x position
 584             if(ch>=2)correctFactor[kBending] = (zMeanChamber[2] - zMeanChamber[0]) / (zRealMatch[2] - zRealMatch[0]);// corrects y position
 585
 586             // searching track intersection with chambers (first approximation)
 587             Float_t deltaZ = zMeanChamber[currCh] - zMeanChamber[0];
 588             trackIntersectCh[currCh][0] = zMeanChamber[currCh] * TMath::Tan(thetaX) * correctFactor[kNonBending];// x position (info from non-bending plane)
 589             trackIntersectCh[currCh][1] = y11 + deltaZ * TMath::Tan(thetaY) * correctFactor[kBending];// y position (info from bending plane)
 590             Int_t detElemIdFromTrack = DetElemIdFromPos(trackIntersectCh[currCh][0], trackIntersectCh[currCh][1], 11+currCh, 0);
 591             if(detElemIdFromTrack<0) {
 592                 if(fDebugLevel>1) printf("Warning: trigger track outside trigger chamber\n");
 593                 continue;
 594             }
 595
 596             triggeredDigits[1][currCh] = MatchingPad(digitStore, detElemIdFromTrack, trackIntersectCh[currCh], matchPad, nboard, zRealMatch, y11);
 597
 598             // if MatchingPad = -500 => too many digits matching pad =>
 599             //                       => Event not clear => Reject track
 600             if(triggeredDigits[1][currCh]<-100){
 601                 isClearEvent = kFALSE;
 602                 if(fDebugLevel>=1) printf("Warning: track = %p (%i) matches many pads. Rejected!\n",recTrigTrack, detElemIdFromTrack);
 603                 break;
 604             }
 605
 606             // deciding if digit matches track
 607             Bool_t isDiffLocBoard = kFALSE;
 608             if(fReproduceTrigResponse && ch>2){
 609                 for(Int_t cath=0; cath<fgkNcathodes; cath++){
 610                     if(boardThatTriggered[currCh][cath][0]>=0){
 611                         if(boardThatTriggered[currCh][cath][0]!=boardThatTriggered[currCh-1][cath][0]) isDiffLocBoard = kTRUE;
 612                     }
 613                 }
 614             }
 615
 616             if(isDiffLocBoard && fDebugLevel>=1)printf("\tDifferent local board\n");
 617
 618             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 619                 match[currCh][cath] = (matchPad[cath] && !isDiffLocBoard);
 620                 if(!match[currCh][cath]) continue;
 621                 digitPerTrack[cath]++;
 622                 trigScheme[currCh][cath]++;
 623                 slatThatTriggered[currCh][cath] = detElemIdFromTrack;
 624                 for(Int_t loc=0; loc<4; loc++){
 625                     boardThatTriggered[currCh][cath][loc] = nboard[cath][loc];
 626                 }
 627             }
 628         } // end chamber loop
 629
 630         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 631             if(digitPerTrack[cath]<3)isClearEvent = kFALSE;
 632             if(fDebugLevel>=1 && !isClearEvent)printf("Warning: found %i digits for trigger track cathode %i.\nRejecting event\n", digitPerTrack[cath],cath);
 633         }
 634
 635         if(!isClearEvent && !fReproduceTrigResponse) continue;
 636
 637         Int_t commonDigits = 0;
 638         for(Int_t ch=0; ch<fgkNchambers; ch++){
 639             if(triggeredDigits[1][ch]==triggeredDigits[0][ch]) commonDigits++; // Compare with previous track
 640             triggeredDigits[0][ch] = triggeredDigits[1][ch]; // Store this track parameters for comparison with next one
 641         }
 642         if(commonDigits>=2){
 643             doubleCountTrack=kTRUE;
 644         }
 645
 646         if(!doubleCountTrack || fReproduceTrigResponse){
 647             for(Int_t cath=0; cath<fgkNcathodes; cath++){
 648                 Int_t is44 = 1;
 649                 Bool_t goodForSlatEff = kTRUE;
 650                 Bool_t goodForBoardEff = kTRUE;
 651                 Int_t ineffSlat = -1;
 652                 Int_t ineffDetElId = -1;
 653                 Int_t firstSlat = slatThatTriggered[0][cath]%100;
 654                 if(firstSlat<0) firstSlat=slatThatTriggered[1][cath]%100;
 655                 Int_t firstBoard = boardThatTriggered[0][kBending][0];
 656                 if(firstBoard<0) firstBoard=boardThatTriggered[1][kBending][0];
 657                 for(Int_t ch=0; ch<fgkNchambers; ch++){
 658                     Bool_t isCurrChIneff = kFALSE;
 659                     is44 *= trigScheme[ch][cath];
 660                     Int_t currSlat = slatThatTriggered[ch][cath]%100;
 661                     if(currSlat<0){
 662                         ineffDetElId = DetElemIdFromPos(trackIntersectCh[ch][0], trackIntersectCh[ch][1], 11+ch, cath);
 663                         currSlat = ineffDetElId%100;
 664                         ineffSlat = currSlat;
 665                         isCurrChIneff = kTRUE;
 666                     }
 667                     if(currSlat!=firstSlat)goodForSlatEff=kFALSE;
 668                     Bool_t atLeastOneLoc=kFALSE;
 669                     if(isCurrChIneff) LocalBoardFromPos(trackIntersectCh[ch][0], trackIntersectCh[ch][1], ineffDetElId, cath, ineffBoard);
 670                     for(Int_t loc=0; loc<4; loc++){
 671                         Int_t currBoard = boardThatTriggered[ch][cath][loc];
 672                         if(isCurrChIneff) currBoard = ineffBoard[loc];
 673                         if(currBoard==firstBoard){
 674                             atLeastOneLoc=kTRUE;
 675                             break;
 676                         }
 677                     }
 678                     if(!atLeastOneLoc)goodForBoardEff=kFALSE;
 679                 } // end chamber loop
 680
 681                 // Trigger 4/4
 682                 if(is44==1){
 683                     fTrigger44[cath]++;
 684                     if(fDebugLevel>=1)printf("Trigger44[%i] = %i\n",cath,fTrigger44[cath]);
 685                     if(goodForSlatEff){
 686                         for(Int_t ch=0; ch<fgkNchambers; ch++){
 687                             fHitPerSlat[ch][cath][firstSlat]++;
 688                             if(fDebugLevel>=1)printf("Slat that triggered = %i\n",slatThatTriggered[ch][cath]);
 689                             if(goodForBoardEff && firstBoard>0){
 690                                 fHitPerBoard[ch][cath][firstBoard-1]++;
 691                                 if(fDebugLevel>=1)printf("Board that triggered = %i\n",firstBoard);
 692                             }
 693                             else if(fDebugLevel>=1) printf("Track = %p: Particle crossed different boards: rejected!\n",recTrigTrack);
 694                         }
 695                     }
 696                     else if(fDebugLevel>=1) printf("Track = %p: Particle crossed different slats: rejected!\n",recTrigTrack);
 697                     //cout<<"fTrigger44["<<cath<<"] = "<<fTrigger44[cath]<<"\tfHitPerSlat["<<0<<"]["<<cath<<"]["<<firstSlat<<"] = "<<fHitPerSlat[0][cath][firstSlat]<<"\tfHitPerBoard["<<0<<"]["<<cath<<"]["<<firstBoard-1<<"] = "<<fHitPerBoard[0][cath][firstBoard-1]<<endl; //REMEMBER TO CUT
 698                 }
 699
 700                 // Trigger 3/4
 701                 if(ineffDetElId>0){
 702                     Int_t ineffCh = ineffDetElId/100-11;
 703                     fTrigger34[ineffCh][cath]++;
 704                     if(fDebugLevel>=1) printf("Trigger34[%i][%i] = %i\n",ineffCh,cath,fTrigger34[ineffCh][cath]);
 705                     if(goodForSlatEff){
 706                         if(fDebugLevel>=1) printf("Slat non efficient = %i\n",ineffDetElId);
 707                         fInefficientSlat[ineffCh][cath][ineffSlat]++;
 708
 709                         if(goodForBoardEff && firstBoard>0){
 710                             if(fDebugLevel>=1) printf("Board non efficient = %i\n",firstBoard);
 711                             fInefficientBoard[ineffCh][cath][firstBoard-1]++;
 712                         }
 713                         else if(fDebugLevel>=1) printf("Track = %p: Particle crossed different boards: rejected!\n",recTrigTrack);
 714                     }
 715                     else if(fDebugLevel>=1) printf("Track = %p: Particle crossed different slats: rejected!\n",recTrigTrack);
 716                     //cout<<"fTrigger34["<<ineffCh<<"]["<<cath<<"] = "<<fTrigger34[ineffCh][cath]<<"\tfInefficientSlat["<<ineffCh<<"]["<<cath<<"]["<<ineffSlat<<"] = "<<fInefficientSlat[ineffCh][cath][ineffSlat]<<"\tfInefficientBoard["<<ineffCh<<"]["<<cath<<"]["<<firstBoard-1<<"] = "<<fInefficientBoard[ineffCh][cath][firstBoard-1]<<endl; //REMEMBER TO CUT
 717                 }
 718             } // end loop on cathodes
 719         }
 720         else if(doubleCountTrack){
 721             if(fDebugLevel>=1)
 722                 printf("\n\tTrack = %p: \nDouble Count Track: Track rejected!\n",recTrigTrack);
 723         }
 724     } // end trigger tracks loop
 725
 726     if(fPrintInfo) InfoDigit(digitStore);
 727 }
 728
 729 //_____________________________________________________________________________
 730 void AliMUONTriggerChamberEff::WriteEfficiencyMap(const char* outputDir)
 731 {
 732     //
 733     /// Writes information on calculated efficiency.
 734     /// It writes: triggerChamberEff.root file containing efficiency histograms.
 735     //
 736
 737     char *cathCode[fgkNcathodes] = {"bendPlane", "nonBendPlane"};
 738
 739     char outFileName[100];
 740
 741     sprintf(outFileName, "%s/MUON.TriggerEfficiencyMap.root",outputDir);
 742     TFile *outputHistoFile = new TFile(outFileName,"RECREATE");
 743     TDirectory *dir = gDirectory;
 744
 745     enum {kSlatIn11, kSlatIn12, kSlatIn13, kSlatIn14, kChamberEff};
 746     char *yAxisTitle = "trigger efficiency (a.u.)";
 747     char *xAxisTitle = "chamber";
 748
 749     TH1F *histo[fgkNcathodes][fgkNchambers+1];
 750     TH1F *histoBoard[fgkNcathodes][fgkNchambers];
 751
 752     // ADDED for check
 753     enum {allChEff, chNonEff, numOfHistoTypes};
 754     char *histoTypeName[numOfHistoTypes] = {"CountInCh", "NonCountInCh"};
 755     char *histoTypeTitle[numOfHistoTypes] = {"counted", "non counted"};
 756     TH1F *histoCheckSlat[fgkNcathodes][fgkNchambers][numOfHistoTypes];
 757     TH1F *histoCheckBoard[fgkNcathodes][fgkNchambers][numOfHistoTypes];
 758     // end ADDED for check
 759
 760     char histoName[40];
 761     char histoTitle[90];
 762
 763     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 764         for(Int_t ch=0; ch<fgkNchambers+1; ch++){
 765             if(ch==kChamberEff){
 766                 sprintf(histoName, "%sChamberEff", cathCode[cath]);
 767                 sprintf(histoTitle, "Chamber efficiency %s", cathCode[cath]);
 768                 histo[cath][ch] = new TH1F(histoName, histoTitle, fgkNchambers, 11-0.5, 15-0.5);
 769                 histo[cath][ch]->SetXTitle(xAxisTitle);
 770                 histo[cath][ch]->SetYTitle(yAxisTitle);
 771                 histo[cath][ch]->GetXaxis()->SetNdivisions(fgkNchambers);
 772             }
 773             else {
 774                 sprintf(histoName, "%sSlatEffChamber%i", cathCode[cath], 11+ch);
 775                 sprintf(histoTitle, "Chamber %i: slat efficiency %s", 11+ch, cathCode[cath]);
 776                 histo[cath][ch] = new TH1F(histoName, histoTitle, fgkNslats, 0-0.5, fgkNslats-0.5);
 777                 histo[cath][ch]->SetXTitle("slat");
 778                 histo[cath][ch]->SetYTitle(yAxisTitle);
 779                 histo[cath][ch]->GetXaxis()->SetNdivisions(fgkNslats);
 780
 781                 sprintf(histoName, "%sBoardEffChamber%i", cathCode[cath], 11+ch);
 782                 sprintf(histoTitle, "Chamber %i: board efficiency %s", 11+ch, cathCode[cath]);
 783                 histoBoard[cath][ch] = new TH1F(histoName, histoTitle, fgkNboards, 1-0.5, fgkNboards+1.-0.5);
 784                 histoBoard[cath][ch]->SetXTitle("boards");
 785                 histoBoard[cath][ch]->SetYTitle(yAxisTitle);
 786                 histoBoard[cath][ch]->GetXaxis()->SetNdivisions(fgkNboards);
 787
 788                 // ADDED for check
 789                 for(Int_t hType=0; hType<numOfHistoTypes; hType++){
 790                     sprintf(histoName, "%sSlat%s%i", cathCode[cath], histoTypeName[hType], 11+ch);
 791                     sprintf(histoTitle, "Chamber %i: slat %s %s", 11+ch, histoTypeTitle[hType], cathCode[cath]);
 792                     histoCheckSlat[cath][ch][hType] = new TH1F(histoName, histoTitle, fgkNslats, 0-0.5, fgkNslats-0.5);
 793                     histoCheckSlat[cath][ch][hType]->SetXTitle("slat");
 794                     histoCheckSlat[cath][ch][hType]->SetYTitle(yAxisTitle);
 795                     histoCheckSlat[cath][ch][hType]->GetXaxis()->SetNdivisions(fgkNslats);
 796
 797                     sprintf(histoName, "%sBoard%s%i", cathCode[cath], histoTypeName[hType], 11+ch);
 798                     sprintf(histoTitle, "Chamber %i: board %s %s", 11+ch, histoTypeTitle[hType], cathCode[cath]);
 799                     histoCheckBoard[cath][ch][hType] = new TH1F(histoName, histoTitle, fgkNboards, 1-0.5, fgkNboards+1.-0.5);
 800                     histoCheckBoard[cath][ch][hType]->SetXTitle("boards");
 801                     histoCheckBoard[cath][ch][hType]->SetYTitle(yAxisTitle);
 802                     histoCheckBoard[cath][ch][hType]->GetXaxis()->SetNdivisions(fgkNboards);
 803                 }
 804                 // end ADDED for check
 805             }
 806         }
 807     }
 808
 809     Float_t efficiency, efficiencyError;
 810     Int_t bin;
 811
 812     for(Int_t cath=0; cath<fgkNcathodes; cath++){
 813         for(Int_t ch=0; ch<fgkNchambers; ch++){
 814             for(Int_t slat=0; slat<fgkNslats; slat++){
 815                 CalculateEfficiency(fHitPerSlat[ch][cath][slat], fHitPerSlat[ch][cath][slat]+fInefficientSlat[ch][cath][slat], efficiency, efficiencyError, kFALSE);
 816                 bin = histo[cath][ch]->FindBin(slat);
 817                 histo[cath][ch]->SetBinContent(bin, efficiency);
 818                 histo[cath][ch]->SetBinError(bin, efficiencyError);
 819
 820                 // ADDED for check
 821                 histoCheckSlat[cath][ch][allChEff]->SetBinContent(bin, fHitPerSlat[ch][cath][slat]);
 822                 histoCheckSlat[cath][ch][chNonEff]->SetBinContent(bin, fInefficientSlat[ch][cath][slat]);
 823             }
 824             CalculateEfficiency(fTrigger44[cath], fTrigger34[ch][cath]+fTrigger44[cath], efficiency, efficiencyError, kFALSE);
 825             bin = histo[cath][ch]->FindBin(11+ch);
 826             histo[cath][kChamberEff]->SetBinContent(bin, efficiency);
 827             histo[cath][kChamberEff]->SetBinError(bin, efficiencyError);
 828
 829             for(Int_t board=0; board<fgkNboards; board++){
 830                 CalculateEfficiency(fHitPerBoard[ch][cath][board], fHitPerBoard[ch][cath][board]+fInefficientBoard[ch][cath][board], efficiency, efficiencyError, kFALSE);
 831                 bin = histoBoard[cath][ch]->FindBin(board+1);
 832                 histoBoard[cath][ch]->SetBinContent(bin, efficiency);
 833                 histoBoard[cath][ch]->SetBinError(bin, efficiencyError);
 834
 835                 // ADDED for check
 836                 histoCheckBoard[cath][ch][allChEff]->SetBinContent(bin, fHitPerBoard[ch][cath][board]);
 837                 histoCheckBoard[cath][ch][chNonEff]->SetBinContent(bin, fInefficientBoard[ch][cath][board]);
 838             }
 839         }
 840     }
 841
 842 // write all histos
 843     outputHistoFile->cd();
 844     dir->GetList()->Write();
 845     outputHistoFile->Close();
 846 }
 847
 848
 849 //_____________________________________________________________________________
 850 void AliMUONTriggerChamberEff::WriteEfficiencyMapTxt(const char* outputDir)
 851 {
 852     //
 853     /// Writes the calculated efficiency in the text file efficiencyCells.dat
 854     ///
 855     /// The file can be further put in $ALICE_ROOT/MUON/data
 856     /// and used to run simulations with measured trigger chamber efficiencies.
 857     //
 858
 859     Int_t effOutWidth=4;
 860
 861     Float_t efficiency, efficiencyError;
 862
 863     Int_t aCapo[] = {16, 38, 60, 76, 92, 108, 117, 133, 155, 177, 193, 209, 225, 234};
 864
 865     char filename[70];
 866     sprintf(filename, "%s/efficiencyCells.dat", outputDir);
 867     ofstream outFile(filename);
 868     outFile << "localBoards" << endl;
 869     for(Int_t ch=0; ch<fgkNchambers; ch++){
 870         //Print information
 871         outFile << "\n\ndetElemId:\t" << 11+ch;
 872         outFile << "00";
 873         for(Int_t cath=0; cath<fgkNcathodes; cath++){
 874             outFile << "\n cathode:\t" << cath << endl;
 875             Int_t currLine=0;
 876             for(Int_t board=0; board<fgkNboards; board++){
 877
 878                 if(board==aCapo[currLine]){
 879                     outFile << endl;
 880                     currLine++;
 881                 }
 882                 CalculateEfficiency(fHitPerBoard[ch][cath][board], fHitPerBoard[ch][cath][board]+fInefficientBoard[ch][cath][board], efficiency, efficiencyError, kFALSE);
 883                 outFile << " " << setw(effOutWidth) << efficiency;
 884             }// loop on boards
 885             outFile << endl;
 886         }// loop on cathodes
 887     }// loop on chambers
 888 }
 889
 890
 891 //_____________________________________________________________________________
 892 Bool_t AliMUONTriggerChamberEff::IsCleanTrack(AliMUONTriggerTrack *triggerTrack,
 893                                               const AliMUONVTrackStore& trackStore)
 894 {
 895     //
 896     /// Try to match track from tracking system with trigger track
 897     //
 898     const Double_t kDistSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY
 899     const Double_t kMaxChi2MatchTrigger = 16.0;
 900
 901     AliMUONTrackParam trackParam;
 902
 903     Double_t distTriggerTrack[3];
 904     Double_t xTrack, yTrack, ySlopeTrack, chi2;
 905
 906     AliMUONTrack* track;
 907     TIter next(trackStore.CreateIterator());
 908
 909     while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
 910     {
 911         trackParam = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->Last()));
 912         AliMUONTrackExtrap::ExtrapToZ(&trackParam, AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber
 913
 914         xTrack = trackParam.GetNonBendingCoor();
 915         yTrack = trackParam.GetBendingCoor();
 916         ySlopeTrack = trackParam.GetBendingSlope();
 917
 918         distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/kDistSigma[0];
 919         distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/kDistSigma[1];
 920         distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/kDistSigma[2];
 921         chi2 = 0.;
 922         for (Int_t iVar = 0; iVar < 3; iVar++) chi2 += distTriggerTrack[iVar]*distTriggerTrack[iVar];
 923         chi2 /= 3.; // Normalized Chi2: 3 degrees of freedom (X,Y,slopeY)
 924         if (chi2 < kMaxChi2MatchTrigger) return kTRUE;
 925     }
 926
 927     return kFALSE;
 928 }
 929
 930
 931 //_____________________________________________________________________________
 932 void AliMUONTriggerChamberEff::SaveInESDFile()
 933 {
 934     //
 935     /// Store AliMUONTriggerChamberEff in esd file
 936     //
 937     TDirectory *dir = gDirectory;
 938     TFile *logFile = 0x0;
 939     TSeqCollection *list = gROOT->GetListOfFiles();
 940     Int_t n = list->GetEntries();
 941     for(Int_t i=0; i<n; i++) {
 942         logFile = (TFile*)list->At(i);
 943         if (strstr(logFile->GetName(), "AliESDs.root")) break;
 944     }
 945     if(logFile){
 946         TTree *esdTree = (TTree*)logFile->Get("esdTree");
 947         if(esdTree){
 948             if(!esdTree->GetUserInfo()->FindObject("AliMUONTriggerChamberEff")){
 949                 AliInfo(Form("Adding AliMUONTrigChamberEff in %s",logFile->GetName()));
 950                 esdTree->GetUserInfo()->Add(this->Clone());
 951                 esdTree->Write("",TObject::kOverwrite);
 952             }
 953         }
 954     }
 955     dir->cd();
 956 }