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