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