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