small fix to remove TF1Helper which sneaked in again, inadvertently...
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliCheckMuonDetEltResponse.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 //Class to check the response of the detection elements of the  MUON tracking chambers 
17 //in function of the position in the detection element.
18 //Author:  Nicolas LE BRIS - SUBATECH Nantes
19 //Modified by Matthieu LENHARDT - SUBATECH Nantes
20
21 //PWG3/muondep:
22 #include "AliCheckMuonDetEltResponse.h"
23
24 //include STEER:
25 #include "AliESDEvent.h"
26 #include "AliESDMuonTrack.h"
27
28 //include MUON:
29 #include "AliMUONTrack.h"
30 #include "AliMUONTrackParam.h"
31 #include "AliMUONTrackExtrap.h"
32 #include "AliMUONVCluster.h"
33 #include "AliMUONConstants.h"
34 #include "AliMUONESDInterface.h"
35 #include "AliMUONGeometryTransformer.h"
36
37 //include MUON/mapping:
38 #include "mapping/AliMpDEManager.h"
39 #include "mapping/AliMpSegmentation.h"
40 #include "mapping/AliMpSlatSegmentation.h"
41 #include "mapping/AliMpSectorSegmentation.h"
42 #include "mapping/AliMpPad.h"
43
44 //include ROOT:
45 #include <TH2F.h>
46 #include <TList.h>
47 #include <TClonesArray.h>
48
49 /// \cond CLASSIMP
50 ClassImp(AliCheckMuonDetEltResponse)
51 /// \endcond
52
53 const Int_t AliCheckMuonDetEltResponse::fgkNCh                   = AliMUONConstants::NTrackingCh();
54 const Int_t AliCheckMuonDetEltResponse::fgkNSt                   = AliMUONConstants::NTrackingSt();
55 const Int_t AliCheckMuonDetEltResponse::fgkNDE                   = 156;
56 const Int_t AliCheckMuonDetEltResponse::fgkNbrOfDetectionElt[10] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26};
57 const Int_t AliCheckMuonDetEltResponse::fgkFirstDetectionElt[10] = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000};
58 const Int_t AliCheckMuonDetEltResponse::fgkOffset                = 100;
59
60 //_____________________________________________________________________________
61 AliCheckMuonDetEltResponse::AliCheckMuonDetEltResponse() 
62 : TObject(),
63   fkTransformer(0x0),
64   fESD(0x0),
65   fTracksTotalNbr(0x0),
66   fNbrUsableTracks(0),
67   fTrackParams(0x0),
68   fTrackParam(0x0),
69   fCluster(0x0),
70   fDetEltTDHistList(0),
71   fDetEltTTHistList(0),
72   fChamberTDHistList(0),
73   fChamberTTHistList(0)
74 {
75 /// Default constructor
76
77     fNbrUsableTracks = 0;
78
79     for (Int_t iCluster = 0; iCluster<fgkNCh; ++iCluster)
80       fNbrClustersCh[iCluster] = 0;
81
82     for (Int_t i=0; i<fgkNCh; ++i)
83       fTrackFilter[i] = 0;
84 }
85
86 //_____________________________________________________________________________
87 AliCheckMuonDetEltResponse::AliCheckMuonDetEltResponse(const AliCheckMuonDetEltResponse& src) 
88 : TObject(src),
89   fkTransformer(0x0),
90   fESD(0x0),
91   fTracksTotalNbr(0x0),
92   fNbrUsableTracks(0),
93   fTrackParams(0x0),
94   fTrackParam(0x0),
95   fCluster(0x0),
96   fDetEltTDHistList(0),
97   fDetEltTTHistList(0),
98   fChamberTDHistList(0),
99   fChamberTTHistList(0)
100 {
101  src.Copy(*this);
102 }
103 //_____________________________________________________________________________
104 AliCheckMuonDetEltResponse& AliCheckMuonDetEltResponse::operator=(const AliCheckMuonDetEltResponse& src) 
105 {
106   /// assignement operator
107   if ( this != &src ) 
108   {
109     src.Copy(*this);
110   }
111   return *this;
112 }
113
114 //_____________________________________________________________________________
115 AliCheckMuonDetEltResponse::AliCheckMuonDetEltResponse(const AliMUONGeometryTransformer* transformer,
116                                                        AliESDEvent* esd,
117                                                        TList* detEltTDHistList,
118                                                        TList* detEltTTHistList,
119                                                        TList* chamberTDHistList,
120                                                        TList* chamberTTHistList) 
121 : TObject(),
122   fkTransformer(transformer),
123   fESD(esd),
124   fTracksTotalNbr(0),
125   fNbrUsableTracks(0),
126   fTrackParams(0x0),
127   fTrackParam(0),
128   fCluster(0),
129   fDetEltTDHistList(detEltTDHistList),
130   fDetEltTTHistList(detEltTTHistList),
131   fChamberTDHistList(chamberTDHistList),
132   fChamberTTHistList(chamberTTHistList)
133 {
134 /// Constructor
135
136     fNbrUsableTracks = 0;
137
138     for (Int_t iCluster = 0; iCluster<fgkNCh; ++iCluster)
139       fNbrClustersCh[iCluster] = 0;
140     
141     for (Int_t i=0; i<fgkNCh; ++i)
142       fTrackFilter[i] = 0;
143 }
144
145
146 //_____________________________________________________________________________
147 AliCheckMuonDetEltResponse::~AliCheckMuonDetEltResponse()
148
149 {
150 /// Destructor
151     delete fTrackParams;
152 }
153
154
155
156 //_____________________________________________________________________________
157 void AliCheckMuonDetEltResponse::CheckDetEltResponse()
158 {
159 //
160 //Cataloging positions (X,Y) of the clusters detected in the detection elements
161 //(fDetEltTDHistList), and positions of crossing points between all the
162 //tracks and the detection elements (fDetEltTTHistList).
163 //Efficiency = 100 * fDetEltTDHistList / fDetEltTTHistList.
164
165 //Loop on tracks
166 //--------------
167   TrackLoop();
168 }
169
170
171
172 //_____________________________________________________________________________
173 void AliCheckMuonDetEltResponse::TrackLoop()
174 {
175 // Check if the track is kept
176   AliESDMuonTrack* esdTrack;
177   AliMUONTrack track;
178   Int_t nTracks, iTrack;
179
180   nTracks = (Int_t)fESD -> GetNumberOfMuonTracks();
181   fTrackParams = new TClonesArray();
182   ///Begininig of the loop:
183   //if (fESD->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL"))
184   {
185     for (iTrack = 0; iTrack < nTracks; iTrack++)
186       {
187         esdTrack   = fESD -> GetMuonTrack(iTrack);
188         
189         if(esdTrack->ContainTrackerData() && esdTrack->GetMatchTrigger() > 0)
190           {
191             AliMUONESDInterface::ESDToMUON(*esdTrack, track);
192             fTrackParams = track.GetTrackParamAtCluster();
193             TrackParamLoop(); //!<Loop on trackParam.
194             fNbrUsableTracks += 1;
195           }
196       }
197   }
198 }
199
200
201
202 //_____________________________________________________________________________
203 void AliCheckMuonDetEltResponse::TrackParamLoop()
204 {
205 // Loop on all the track params and fill the histos
206   Int_t nTrackParams = (Int_t) fTrackParams->GetEntriesFast();  //!<Number of trackParams in the track.
207   Int_t iTrackParam = 0;                                        //!<Number of the trackParam of the track.
208   Int_t oldChamber = -1, newChamber = 0; //!<To check if there is 0, 1 or 2 (overlap cases) clusters in the same chamber for a track.                                      
209   Int_t detElt;                          //!<Detection element Id.
210   
211   for (Int_t ch = 0; ch < fgkNCh; ++ch)
212     fTrackFilter[ch] = 0;
213
214   Double_t posXL, posYL, posZL;          //!<Local positions.
215   Double_t posXG, posYG, posZG;          //!<Global. positions.
216   Int_t chamberResponse [10] = {0};      //!<1 if the chamber has responded; 0 if not
217   
218   for (iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
219     { 
220       fTrackParam = (AliMUONTrackParam*) fTrackParams->At(iTrackParam);
221       fCluster    = (AliMUONVCluster*  ) fTrackParam ->GetClusterPtr();    
222       fTrackFilter   [fCluster->GetChamberId()] = 1;
223       chamberResponse[fCluster->GetChamberId()] = 1;
224     }
225
226   for (Int_t station = 0; station < fgkNSt-1; ++station)
227     {
228       Int_t filter;                                                       //<!
229       Int_t ch1, ch2, ch3, ch4;                                           //<!
230       ch1 = 2*station;                                                    //<!
231       ch2 = 2*station + 1;                                                //<!
232       ch3 = 2*station + 2;                                                //<!
233       ch4 = 2*station + 3;                                                //<!
234                                                                           //<!For the efficiency calculation the tracks
235       if (station < 3 )                                                   //<!reconstructed must have responded to the
236         {                                                                 //<!criteria of the tracking. 
237           filter            = fTrackFilter[ch1];                          //<!And that's why the tracks usable for the 
238           fTrackFilter[ch1] = fTrackFilter[ch2];                          //<!intrinsic efficiency calculation are
239           fTrackFilter[ch2] = filter;                                     //<!the tracks which have one or two clusters
240         }                                                                 //<!in each station. So the case where a track
241                                                                           //<!hasn't a cluster in a station is not
242       else                                                                //<!taking into account.
243         {                                                                 //<!This part solves the problem. See the ALICE 
244           if (chamberResponse[ch3]*chamberResponse[ch4] != 0)             //<!note of Diego STOCCO on the trigger efficiency
245             {                                                             //<!
246               filter            = fTrackFilter[ch1];                      //<!
247               fTrackFilter[ch1] = fTrackFilter[ch2];                      //<!
248               fTrackFilter[ch2] = filter;                                 //<!
249             }                                                             //<!
250           else                                                            //<!
251             {                                                             //<!
252               fTrackFilter[ch1] = 0;                                      //<!
253               fTrackFilter[ch2] = 0;                                      //<!
254             }                                                             //<!
255
256           if (chamberResponse[ch1]*chamberResponse[ch2] != 0)
257             {
258               filter            = fTrackFilter[ch3];
259               fTrackFilter[ch3] = fTrackFilter[ch4];
260               fTrackFilter[ch4] = filter;
261             }
262           else
263             {
264               fTrackFilter[ch3] = 0;
265               fTrackFilter[ch4] = 0;
266             }
267         }
268     }
269   
270
271   ///Begining of the loop:
272   for (iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
273     {
274       fTrackParam = (AliMUONTrackParam*) fTrackParams->At(iTrackParam);
275       fCluster    = (AliMUONVCluster*  ) fTrackParam ->GetClusterPtr(); 
276       
277       newChamber  = fCluster->GetChamberId();
278       detElt      = fCluster->GetDetElemId();
279
280       ///Global and local positions calculation:
281       posXG = fTrackParam->GetNonBendingCoor(); 
282       posYG = fTrackParam->GetBendingCoor(); 
283       posZG = fTrackParam->GetZ(); 
284       
285       fkTransformer->Global2Local(detElt, posXG, posYG, posZG, posXL, posYL, posZL);  //!<Transfomation from global to local positions.
286       
287       ///Filling histograms of the cluster positions on the detection element of the TRACKS DETECTED (TD):
288       FillTDHistos(newChamber, detElt, posXL, posYL);
289     
290       ///Filling histograms of the cluster positions on the detection element of ALL THE TRACKS (TT):
291       FillTTHistos(newChamber, detElt, posXL, posYL);
292
293       if (newChamber != oldChamber) 
294         {
295           if (newChamber > oldChamber + 1)                                 //!<Check if it doesn't miss a chamber.
296             {
297               Int_t nbrMissChamber = newChamber - (oldChamber + 1);
298               FindAndFillMissedDetElt(fTrackParam, oldChamber+1, nbrMissChamber); //!<Calculation of the parameters of the missing cluster(s).
299             }
300             
301           if ( iTrackParam == nTrackParams - 1 && newChamber != fgkNCh-1)           //!<Check if the last chamber, chamber 9 (from 0 to 9) has responded.
302             FindAndFillMissedDetElt(fTrackParam, fgkNCh-1, 1);                      //!<Calculation of the parameters of the missing cluster(s) in the last chamber.
303             
304         }
305       oldChamber = newChamber; 
306     } 
307 }
308
309
310
311 //_____________________________________________________________________________
312 void AliCheckMuonDetEltResponse::FillTDHistos(Int_t chamber,
313                                               Int_t detElt,
314                                               Double_t posXL,
315                                               Double_t posYL)
316 {
317 // Fill the histo for tracks detected
318   if(fTrackFilter[chamber]== 1)
319     {
320       Int_t iDet = 0; //!<Position of the detection element in the histograms' list.
321       iDet = FromDetElt2iDet(chamber, detElt);
322       ((TH2F*) fDetEltTDHistList->At(iDet))->Fill(posXL, posYL);
323
324       Int_t detEltLocalId = 0;  //!<Id of the detection element in the station
325       detEltLocalId =  FromDetElt2LocalId(chamber, detElt);
326       ((TH1F*) fChamberTDHistList->At(chamber))->Fill(detEltLocalId);
327       ((TH1F*) fChamberTDHistList->At(fgkNCh))->Fill(chamber + 1);
328    }
329 }
330
331
332
333
334 //_____________________________________________________________________________
335 void AliCheckMuonDetEltResponse::FillTTHistos(Int_t chamber,
336                                               Int_t detElt,
337                                               Double_t posXL,
338                                               Double_t posYL)
339 {
340 // Fill the histo for total number of tracks
341   if(fTrackFilter[chamber] == 1)
342     {
343       Int_t iDet = 0; //!<Position of the detection element in the histograms' list.
344       iDet = FromDetElt2iDet(chamber, detElt);
345       ((TH2F*) fDetEltTTHistList->At(iDet)) -> Fill(posXL, posYL);
346      
347       Int_t detEltLocalId = 0;  //!<Id of the detection element in the station
348       detEltLocalId =  FromDetElt2LocalId(chamber, detElt);
349       ((TH1F*) fChamberTTHistList->At(chamber))->Fill(detEltLocalId);
350       ((TH1F*) fChamberTTHistList->At(fgkNCh))->Fill(chamber + 1);
351     }
352 }
353
354
355
356
357
358 //_____________________________________________________________________________
359 Int_t AliCheckMuonDetEltResponse::FromDetElt2iDet(Int_t chamber, 
360                                                   Int_t detElt) const
361 {
362 ///
363 ///Connexion between the detection element X and its position in the list of histograms iX.
364 ///
365
366     Int_t iDet = 0; //!<Position of the detection element (detElt) in the histograms' list.
367
368     if (chamber<4)             iDet = detElt-fgkOffset*(chamber+1)+ 4* chamber      ; 
369     if (chamber>3 && chamber<6) iDet = detElt-fgkOffset*(chamber+1)+18*(chamber-4)+16;
370     if (chamber>5)             iDet = detElt-fgkOffset*(chamber+1)+26*(chamber-6)+52;
371
372     return iDet;    
373 }
374
375
376
377 //_____________________________________________________________________________
378 Int_t AliCheckMuonDetEltResponse::FromDetElt2LocalId(Int_t chamber, 
379                                                      Int_t detElt) const
380 {
381 ///
382 ///Connexion between the detection element X and its number in the station.
383 ///
384
385     Int_t localId = 0; //!<Position of the detection element (detElt) in the histograms' list.
386     localId = detElt - (chamber+1) * 100;
387
388     return localId;    
389 }
390
391
392
393 //_____________________________________________________________________________
394 void AliCheckMuonDetEltResponse::FindAndFillMissedDetElt(AliMUONTrackParam* extrapTrackParam, 
395                                                          Int_t firstMissCh,
396                                                          Int_t nbrMissCh)
397 {
398 ///
399 ///Find which detection elements should have been hit but were missed, 
400 ///and fill the TT histos appropriately
401 ///
402   for (Int_t iCh = 0; iCh < nbrMissCh; ++iCh)
403     {
404       Int_t chamber = firstMissCh + iCh;
405       Int_t nbrOfDetElt =  AliMpDEManager::GetNofDEInChamber(chamber, kTRUE); //!<Number of detection elements in the chamber.
406       
407       Double_t pos1[6] = {0, 0, 0, 0, 0, 0};        //!<First point used to compute the extrapolated point (first 3 for global coordinates, last 3 for local).
408       Double_t pos2[6] = {0, 0, 0, 0, 0, 0};        //!<Second point used to compute the extrapolated point (first 3 for global coordinates, last 3 for local).
409       Double_t posMiss[2] = {0, 0};                 //!<(X, Y) local coordinates of the missing cluster.
410             
411       pos1[2] = AliMUONConstants::DefaultChamberZ(chamber);           //!<Z of point 1, defined by being the Z of the chamber in "perfect" position.
412       AliMUONTrackExtrap::ExtrapToZ(extrapTrackParam, pos1[2]);
413       pos1[0] = extrapTrackParam->GetNonBendingCoor();                //!<X of point 1, extrapolated by following the Track.
414       pos1[1] = extrapTrackParam->GetBendingCoor();                   //!<Y of point 1, extrapolated by following the Track.
415       
416       pos2[2] = AliMUONConstants::DefaultChamberZ(chamber) + AliMUONConstants::DzCh();   //!<Z of point 2, defined by being the Z of the chamber in "perfect" position 
417       AliMUONTrackExtrap::ExtrapToZ(extrapTrackParam, pos2[2]);                           //!< + plus a small shift (the distance between two stations in a same chamber).
418       pos2[0] = extrapTrackParam->GetNonBendingCoor();                                   //!<X of point 2, extrapolated by following the Track.                        
419       pos2[1] = extrapTrackParam->GetBendingCoor();                                      //!<Y of point 2, extrapolated by following the Track.                          
420       
421       
422       
423         for (Int_t iDE = 0; iDE < nbrOfDetElt; iDE++)                    //!<Loop on all the detection element of the chamber
424           {
425             Int_t deId = (chamber + 1)*fgkOffset + iDE;                        //!<detection element Id 
426             
427             fkTransformer->Global2Local(deId, pos1[0], pos1[1], pos1[2], pos1[3], pos1[4], pos1[5]);      //!<convesrion of point 1 and 2 in the local coordinates
428             fkTransformer->Global2Local(deId, pos2[0], pos2[1], pos2[2], pos2[3], pos2[4], pos2[5]);
429             
430             CoordinatesOfMissingCluster(pos1[3], pos1[4], pos1[5], pos2[3], pos2[4], pos2[5], posMiss[0], posMiss[1]);
431
432             Bool_t isMissed = kFALSE;
433             if (chamber < 4)
434               isMissed = CoordinatesInDetEltSt12(deId, posMiss[0], posMiss[1]);
435             else
436               isMissed = CoordinatesInDetEltSt345(deId, posMiss[0], posMiss[1]);
437
438             if (isMissed)
439               FillTTHistos(chamber, deId, posMiss[0], posMiss[1]);
440           }
441     }
442 }
443
444
445
446 //_____________________________________________________________________________
447 void AliCheckMuonDetEltResponse::CoordinatesOfMissingCluster(Double_t x1, Double_t y1, Double_t z1,
448                                                              Double_t x2, Double_t y2, Double_t z2,
449                                                              Double_t& x, Double_t& y) const
450 {
451   //
452   //Compute the coordinates of the missing cluster.
453   //There are defined by the intersection between the straigth line joining two extrapolated points (1 and 2) and the detection element plane.
454   //In the local coordinates, this means Z=0 in the parametric equation of the line.
455   //
456
457   Double_t t = 0;
458   t = - z1 / (z2 - z1);
459   
460   x = t * (x2 - x1) + x1;
461   y = t * (y2 - y1) + y1;
462 }
463
464
465 //_____________________________________________________________________________
466 Bool_t AliCheckMuonDetEltResponse::CoordinatesInDetEltSt345(Int_t DeId, Double_t x, Double_t y)
467 {
468   //
469   //Return kTRUE if the coordinates are in the Detection Element, for station 3, 4 and 5.
470   //This is done by checking if a pad correspond to the (x, y) position.
471   //  
472
473   AliMpPad pad1;
474   AliMpPad pad2;
475
476   AliMpSlatSegmentation *segm1 = new AliMpSlatSegmentation(AliMpSegmentation::Instance(kFALSE)->GetSlat(DeId, AliMp::kCath0));
477   AliMpSlatSegmentation *segm2 = new AliMpSlatSegmentation(AliMpSegmentation::Instance(kFALSE)->GetSlat(DeId, AliMp::kCath1));
478   pad1 = segm1->PadByPosition(x, y, kFALSE);
479   pad2 = segm2->PadByPosition(x, y, kFALSE);
480  
481   if (pad1.IsValid() && pad2.IsValid())
482     return kTRUE;
483   else
484     return kFALSE;
485 }
486
487
488 //_____________________________________________________________________________
489 Bool_t AliCheckMuonDetEltResponse::CoordinatesInDetEltSt12(Int_t DeId, Double_t x, Double_t y)
490 {
491   //Return kTRUE if the coordinates are in the Detection Element, for station 1 and 2.
492   //This is done by checking if a pad correspond to the (x, y) position.
493   
494   AliMpPad pad1;
495   AliMpPad pad2;
496
497   AliMpSectorSegmentation *segm1 = new AliMpSectorSegmentation(AliMpSegmentation::Instance(kFALSE)->GetSector(DeId, AliMp::kCath0));
498   AliMpSectorSegmentation *segm2 = new AliMpSectorSegmentation(AliMpSegmentation::Instance(kFALSE)->GetSector(DeId, AliMp::kCath1));
499   pad1 = segm1->PadByPosition(x, y, kFALSE);
500   pad2 = segm2->PadByPosition(x, y, kFALSE);
501  
502   if (pad1.IsValid() && pad2.IsValid())
503     return kTRUE;
504   else
505     return kFALSE;
506 }