end-of-line normalization
[u/mrichter/AliRoot.git] / ITS / UPGRADE / KMCDetector.h
1 #ifndef DETECTORK_H
2 #define DETECTORK_H
3
4 #include <TNamed.h>
5 #include <TList.h>
6 #include <TGraph.h>
7 #include <Riostream.h>
8 #include <TMatrixD.h>
9 #include <TClonesArray.h>
10 #include <TArrayD.h>
11 #include <TH2.h>
12 #include <TH1.h>
13 #include "AliExternalTrackParam.h"
14 #include "AliLog.h"
15
16 //-------------------------------------------------------------------------
17 // Current support and development: Ruben Shahoyan (Ruben.Shahoyan@cern.ch) 
18 //-------------------------------------------------------------------------
19
20
21 class KMCLayer;
22 class KMCCluster;
23
24 //////////////////////////////////////////////////////////////////////////////////////////////
25 //--------------------------------------------------------------------------------------------
26 class KMCProbe : public AliExternalTrackParam {
27  public:
28   enum {kBitKilled=BIT(14)};
29   enum {kNDOF=5};
30   enum {kY2=0,kZ2=2,kSnp2=5,kTgl2=9,kPtI2=14};
31   enum {kY,kZ,kSnp,kTgl,kPtI};
32   //
33  KMCProbe() : fMass(0.14),fChi2(0),fHits(0),fFakes(0),fNHits(0),fNHitsITS(0),fNHitsITSFake(0),fInnLrCheck(fgNITSLayers) {for (int i=kMaxITSLr;i--;) fClID[i] =-2;
34 }
35   KMCProbe(KMCProbe& src);
36   KMCProbe& operator=(const KMCProbe& src);
37   virtual ~KMCProbe() {}
38   virtual   void   Print(Option_t* option = "") const;
39   virtual   void   Reset() {fMass=0.14; fChi2=0; fHits=fFakes=0;  for (int i=kMaxITSLr;i--;) fClID[i]=-2; AliExternalTrackParam::Reset();}
40   virtual   Bool_t IsSortable()                     const {return kTRUE;}
41   virtual   Int_t  Compare(const TObject* obj)      const;
42   void      ResetCovMat();
43   //
44   void      Kill(Bool_t v=kTRUE)                          {SetBit(kBitKilled,v);}
45   Bool_t    IsKilled()                              const {return TestBit(kBitKilled);}
46   //
47   Bool_t    CorrectForMeanMaterial(const KMCLayer* lr, Bool_t inward=kTRUE);
48   Bool_t    GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir=0) const;
49   Bool_t    PropagateToR(double r, double b, int dir=0);
50   
51   //
52   void      SetMass(double m=0.14)                      {fMass = m;}
53   Double_t  GetMass()                             const {return fMass;}
54   Double_t  GetChi2()                             const {return fChi2;}
55   void      SetChi2(double chi2)                        {fChi2 = chi2;}
56   void      AddChi2(double chi2)                        {fChi2 += chi2;}
57   void      SetInnerLrChecked(Int_t n)                  {if (n<fgNITSLayers) fInnLrCheck = n;}
58   UShort_t  GetNITSHits()                         const {return fNHitsITS;}
59   UShort_t  GetNFakeITSHits()                     const {return fNHitsITSFake;}
60   UShort_t  GetNHits()                            const {return fNHits;}
61   UShort_t  GetInnerLayerChecked()                const {return fInnLrCheck;}
62   UInt_t&   GetHitsPatt()                               {return fHits;}
63   UInt_t&   GetFakesPatt()                              {return fFakes;}
64   Double_t  GetNormChi2(Bool_t penalize=kFALSE)   const;
65   void      AddHit(Int_t lr, double chi2, Int_t clID=-1);
66   void      ResetHit(Int_t lr);
67   Bool_t    IsHit(Int_t lr)                       const {return (lr<fgNITSLayers) ? IsWBit(fHits,lr)  : kFALSE;}
68   Bool_t    IsHitFake(Int_t lr)                   const {return (lr<fgNITSLayers) ? IsWBit(fFakes,lr) : kFALSE;}
69   Bool_t    PropagateToCluster(KMCCluster* cl, Double_t b);
70   // protected: 
71   static void   SetWBit(UInt_t &patt,UInt_t bit)               {patt |= 0x1<<bit;}
72   static void   ResetWBit(UInt_t &patt,UInt_t bit)             {patt &= ~(0x1<<bit);}
73   static Bool_t IsWBit(const UInt_t &patt,const UInt_t bit)    {return patt&(0x1<<bit);}
74   static void   SetNITSLayers(Int_t n)                         {fgNITSLayers = n;}
75   static int    GetNITSLayers()                                {return fgNITSLayers;}
76   //
77   static Double_t GetMissingHitPenalty()                        {return fgMissingHitPenalty;}
78   static void     SetMissingHitPenalty(double p=2.)             {fgMissingHitPenalty = p;}
79   //
80  public:
81   enum {kMaxITSLr=12};
82   Float_t fMass;   // mass
83   Float_t fChi2;   // total chi2
84   UInt_t  fHits;   // pattern on hits (max 32!)
85   UInt_t  fFakes;  // pattern of fakes among hits
86   UShort_t fNHits;    // total hits
87   UShort_t fNHitsITS; // total ITS hits
88   UShort_t fNHitsITSFake; // number of fake ITS hits
89   UShort_t fInnLrCheck;   // lowest active layer where update was checked
90   Int_t    fClID[kMaxITSLr];
91   //
92   static Int_t    fgNITSLayers;
93   static Double_t fgMissingHitPenalty;
94   ClassDef(KMCProbe,2);  
95 };
96
97 //_______________________________________
98 inline Double_t KMCProbe::GetNormChi2(Bool_t penalize) const
99 {
100   // normalized chi2, penilized for missing hits
101   //  if (fNHitsITS<3) return 0;
102   double chi2 = fChi2;
103   if (penalize) {
104     int nMiss = fgNITSLayers - fNHitsITS - fInnLrCheck;
105     chi2 = fChi2 + fgMissingHitPenalty*nMiss;
106   }
107   return chi2/( (fNHitsITS<<1)-kNDOF);
108 }
109
110 //_______________________________________
111 inline void KMCProbe::AddHit(Int_t lr, double chi2, Int_t clID) {
112   // note: lr is active layer ID
113   if (lr<0) return;
114   fNHits++;
115   if (lr<fgNITSLayers) {
116     fChi2 += chi2;
117     SetWBit(fHits,lr); 
118     fNHitsITS++;
119     if (clID>-1) {
120       SetWBit(fFakes,lr);
121       fNHitsITSFake++;
122     }
123     fClID[lr] = clID;
124     //else ResetWBit(fFakes,lr);
125   }
126 }
127
128 //_______________________________________
129 inline void KMCProbe::ResetHit(Int_t lr) {
130   // note: lr is active layer ID
131   if (lr>=fgNITSLayers) return; 
132   if (IsWBit(fHits,lr))  {fNHitsITS--;     ResetWBit(fHits,lr);}
133   if (IsWBit(fFakes,lr)) {fNHitsITSFake--; ResetWBit(fFakes,lr);}
134 }
135
136 //////////////////////////////////////////////////////////////////////////////////////////////
137 //--------------------------------------------------------------------------------------------
138 class KMCCluster : public TObject {
139  public:
140   //
141   enum {kBitKilled=BIT(14)};
142  KMCCluster(Float_t y=0, Float_t z=0, Float_t x=0, Float_t phi=0) : fY(y),fZ(z),fX(x),fPhi(phi) {}
143   KMCCluster(KMCCluster &src);
144   KMCCluster& operator=(const KMCCluster& src);
145   virtual ~KMCCluster() {}
146   void Reset()               {Clear();}
147   //
148   Double_t GetY()    const   {return fY;}
149   Double_t GetX()    const   {return fX;}
150   Double_t GetZ()    const   {return fZ;}
151   Double_t GetPhi()  const   {return fPhi;}
152   //
153   void    Kill(Bool_t v=kTRUE)          {SetBit(kBitKilled,v);}
154   Bool_t  IsKilled()              const {return TestBit(kBitKilled);}
155   Float_t fY; 
156   Float_t fZ; 
157   Float_t fX;
158   Float_t fPhi;
159   void Set(Float_t y, Float_t z, Float_t x, Float_t phi) {fY=y; fZ=z; fX=x; fPhi=phi; ResetBit(kBitKilled);}
160   virtual void Print(Option_t * = 0) const;
161   ClassDef(KMCCluster,1);
162 };
163
164 //////////////////////////////////////////////////////////////////////////////////////////////
165 //--------------------------------------------------------------------------------------------
166 class KMCLayer : public TNamed {
167 public:
168   enum {kTypeNA=-1,kITS,kTPC,kBitVertex=BIT(15)};
169   KMCLayer(char *name);
170   Float_t GetRadius()   const {return fR;}
171   Float_t GetRadL()     const {return fx2X0;}
172   Float_t GetXTimesRho() const {return fXRho;}
173   Float_t GetPhiRes()   const {return fPhiRes;}
174   Float_t GetZRes()     const {return fZRes;}
175   Float_t GetLayerEff() const {return fEff;}
176   Int_t   GetActiveID() const {return fActiveID;}
177   virtual void  Print(Option_t* option = "") const;
178   //
179   Bool_t IsDead()       const {return fIsDead;}
180   Bool_t IsITS()        const {return fType==kITS;}
181   Bool_t IsTPC()        const {return fType==kTPC;}
182   Bool_t IsVertex()     const {return TestBit(kBitVertex);}
183   //
184   Int_t    AddBgCluster(double y,double z,double x,double phi) {int n=GetNBgClusters(); new (fClBg[n]) KMCCluster(y,z,x,phi); return ++n;}
185   KMCCluster* GetBgCluster(Int_t i) {return (KMCCluster*)fClBg[i];}
186   KMCCluster* GetMCCluster()        {return (KMCCluster*)&fClMC;}
187   //
188   KMCProbe*       GetAnProbe()    const {return (KMCProbe*)&fTrCorr;}
189   Int_t           GetNMCTracks()  const {return fTrMC.GetEntries();}
190   KMCProbe*       AddMCTrack(KMCProbe* src=0);
191   KMCProbe*       GetMCTrack(Int_t it) const {return (KMCProbe*)fTrMC[it];}
192   KMCProbe*       GetWinnerMCTrack()  {if (!fTrMC.IsSorted()) fTrMC.Sort(); return fTrMC.GetEntries() ? (KMCProbe*)fTrMC[0]:0;}
193   TClonesArray*         GetMCTracks()        const {return (TClonesArray*)&fTrMC;}
194   void Reset();
195   void ResetMC() { fClBg.Clear(); fTrMC.Clear();}
196   void ResetBgClusters() { fClBg.Clear(); }
197   void ResetMCTracks()   { fTrMC.Clear(); }
198   Int_t GetNBgClusters() const { return fClBg.GetEntries();}
199   static Double_t GetDefEff()   {return fgDefEff;}
200   static void     SetDefEff(double eff=1) {fgDefEff = eff>1. ? 1.: (eff<0? 0:eff);}
201   //
202   //
203   Float_t fR; 
204   Float_t fx2X0;
205   Float_t fXRho;    // x*density
206   Float_t fPhiRes; 
207   Float_t fZRes;   
208   Float_t fEff;
209   Bool_t  fIsDead;
210   Int_t   fType;   // its, tpc etc
211   Int_t   fActiveID;   // active layer id
212   Float_t fSig2EstD;
213   Float_t fSig2EstZ;
214   //
215   KMCCluster   fClCorr;     // ideal cluster
216   KMCCluster   fClMC;       // MC cluster (from MS scattered track)
217   TClonesArray fClBg; // bg clusters for MC
218   //
219   KMCProbe     fTrCorr;   // ideal track
220   TClonesArray fTrMC;      // MC tracks
221   //
222   static Double_t fgDefEff;
223   ClassDef(KMCLayer,1);
224 };
225
226 //////////////////////////////////////////////////////////////////////////////////////////////
227 //--------------------------------------------------------------------------------------------
228 class KMCDetector : public TNamed {
229  public:
230   enum {kUtilHisto=BIT(14)};
231   KMCDetector();
232   KMCDetector(char *name,char *title);
233   virtual ~KMCDetector();
234
235   void AddLayer(char *name, Float_t radius, Float_t radL, Float_t xrho=0, Float_t phiRes=999999, Float_t zRes=999999, Float_t eff=-1);
236   Int_t GetLayerID(Int_t actID) const;
237   void KillLayer(char *name);
238   void SetRadius(char *name, Float_t radius);
239   void SetRadiationLength(char *name, Float_t radL);
240   void SetResolution(char *name, Float_t phiRes=999999, Float_t zRes=999999);
241   void SetLayerEfficiency(char *name, Float_t eff=1.0);
242   void RemoveLayer(char *name);
243
244   Float_t GetRadius(char *name);
245   Float_t GetRadiationLength(char *name);
246   Float_t GetResolution(char *name, Int_t axis=0);
247   Float_t GetLayerEfficiency(char *name);
248
249   void PrintLayout(); 
250   void PlotLayout(Int_t plotDead = kTRUE);
251   
252   void MakeAliceAllNew(Bool_t flagTPC =0,  Bool_t flagMon=1,int setVer=0);
253   void MakeAliceCurrent(Bool_t flagTPC =0, Int_t AlignResiduals = 0);
254   void AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip=1);
255   void RemoveTPC();
256
257   void SetBField(Float_t bfield) {fBFieldG = bfield*10; }
258   Float_t GetBField() const {return fBFieldG/10; }
259   void SetLhcUPCscale(Float_t lhcUPCscale) {fLhcUPCscale = lhcUPCscale; }
260   Float_t GetLhcUPCscale() const { return fLhcUPCscale; }
261   void SetParticleMass(Float_t particleMass) {fParticleMass = particleMass; }
262   Float_t GetParticleMass() const { return fParticleMass; }
263   void SetIntegrationTime(Float_t integrationTime) {fIntegrationTime = integrationTime; }
264   Float_t GetIntegrationTime() const { return fIntegrationTime; }
265   void SetMaxRadiusOfSlowDetectors(Float_t maxRadiusSlowDet) {fMaxRadiusSlowDet =  maxRadiusSlowDet; }
266   Float_t GetMaxRadiusOfSlowDetectors() const { return fMaxRadiusSlowDet; }
267   void SetAvgRapidity(Float_t avgRapidity) {fAvgRapidity = avgRapidity; }
268   Float_t GetAvgRapidity() const { return fAvgRapidity; }
269   void SetConfidenceLevel(Float_t confLevel) {fConfLevel = confLevel; }
270   Float_t GetConfidenceLevel() const { return fConfLevel; }
271   void SetAtLeastCorr(Int_t atLeastCorr ) {fAtLeastCorr = atLeastCorr; }
272   Float_t GetAtLeastCorr() const { return fAtLeastCorr; }
273   void SetAtLeastFake(Int_t atLeastFake ) {fAtLeastFake = atLeastFake; }
274   Float_t GetAtLeastFake() const { return fAtLeastFake; }
275
276   void SetdNdEtaCent(Int_t dNdEtaCent ) {fdNdEtaCent = dNdEtaCent; }
277   Float_t GetdNdEtaCent() const { return fdNdEtaCent; }
278
279   Int_t GetNLayers()          const {return fLayers.GetEntries(); }
280   Int_t GetNActiveLayers()    const {return fNActiveLayers; }
281   Int_t GetNActiveITSLayers() const {return fNActiveITSLayers; }
282
283   Bool_t SolveSingleTrackViaKalman(Double_t mass, Double_t pt, Double_t eta);
284   Bool_t SolveSingleTrackViaKalmanMC(int offset=6);
285   Bool_t SolveSingleTrack(Double_t mass, Double_t pt, Double_t eta, TObjArray* sumArr=0, int nMC=10000,int offset=6);
286   KMCProbe* KalmanSmooth(int actLr, int actMin,int actMax) const;
287   KMCProbe* KalmanSmoothFull(int actLr, int actMin,int actMax) const; //TBD
288   void   EliminateUnrelated();
289   //
290   KMCProbe* PrepareKalmanTrack(double pt, double lambda, double mass, int charge, double phi=0,double x=0,double y=0,double z=0);
291   Bool_t TransportKalmanTrackWithMS(KMCProbe *probTr, int maxLr) const;
292   void   ApplyMS(KMCProbe* trc,  double x2x0) const;
293   Bool_t PropagateToLayer(KMCProbe* trc, KMCLayer* lr, int dir) const;
294   Bool_t UpdateTrack(KMCProbe* trc, KMCLayer* lr, KMCCluster* cl, Bool_t goToCluster=kTRUE) const;
295
296   // Helper functions
297   Double_t ThetaMCS                 ( Double_t mass, Double_t RadLength, Double_t momentum ) const;
298   Double_t ProbGoodHit              ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; 
299   Double_t ProbGoodChiSqHit         ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; 
300   Double_t ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; 
301   Double_t ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; 
302  
303   // Howard W. hit distribution and convolution integral
304   Double_t Dist              ( Double_t Z, Double_t radius ) ;  
305   Double_t HitDensity        ( Double_t radius )   ;
306   Double_t UpcHitDensity     ( Double_t radius )   ;
307   Double_t IntegratedHitDensity  ( Double_t multiplicity, Double_t radius )   ;
308   Double_t OneEventHitDensity    ( Double_t multiplicity, Double_t radius ) const   ;
309   Double_t D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][20] ) const ;
310   
311   TGraph* GetGraphMomentumResolution(Int_t color, Int_t linewidth=1);
312   TGraph* GetGraphPointingResolution(Int_t axis,Int_t color, Int_t linewidth=1);
313   TGraph* GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth=1);
314
315   TGraph* GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth=1);
316
317   TGraph* GetGraphRecoEfficiency(Int_t particle, Int_t color, Int_t linewidth=1); 
318   TGraph* GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth);
319
320   TGraph* GetGraph(Int_t number, Int_t color, Int_t linewidth=1);
321
322   void MakeStandardPlots(Bool_t add =0, Int_t color=1, Int_t linewidth=1,Bool_t onlyPionEff=0);
323
324   // method to extend AliExternalTrackParam functionality
325   Bool_t IsZero(double val, double tol=1e-9) const {return TMath::Abs(val)<tol;}
326   TList *GetLayers()                   const {return (TList*)&fLayers;}
327   KMCLayer* GetLayer(Int_t i)          const {return (KMCLayer*) fLayers.At(i);}
328   KMCLayer* GetActiveLayer(Int_t actID)    const {int pid=GetLayerID(actID); return pid<0 ? 0:GetLayer(pid);}
329   KMCLayer* GetLayer(const char* name) const {return (KMCLayer*) fLayers.FindObject(name);}
330   KMCProbe* GetProbeTrack()       const {return (KMCProbe*)&fProbe;}
331   void   ClassifyLayers();
332   void   ResetMCTracks(int maxLr=-1);                      
333   //
334   Bool_t   GetUseBackground()               const {return fUseBackground;}
335   void     SetUseBackground(Bool_t v=kTRUE)       {fUseBackground = v;}
336   void     CheckTrackProlongations(KMCProbe *probe, KMCLayer* lr, KMCLayer* lrP);
337   void     ResetSearchLimits() {fBgYMin=fBgZMin=1e6; fBgYMax=fBgZMax=-1e6; fNBgLimits=0;}
338   void     UpdateSearchLimits(KMCProbe* probe, KMCLayer* lr);
339   Int_t    GenBgClusters(KMCLayer* lr);
340   Bool_t   NeedToKill(KMCProbe* probe) const;
341   Double_t PropagateBack(KMCProbe* trc);
342   //
343   // definition of reconstructable track
344   void     RequireMaxChi2Cl(double cut=25.)           {fMaxChi2Cl = cut>0 ? cut:9; fMaxChi2ClSQ = TMath::Sqrt(fMaxChi2Cl);}
345   void     RequireMinITSHits(Int_t n=4)               {fMinITSHits = n;}
346   void     RequireMaxNormChi2NDF(double cut=5.)       {fMaxNormChi2NDF = cut>0 ? cut:9;}
347   void     RequirePattern(UInt_t *patt, int groups);
348   //
349   Double_t GetMaxChi2Cl()                      const {return fMaxChi2Cl;}
350   Double_t GetMaxNormChi2NDFusterKMC()              const {return fMaxNormChi2NDF;}
351   Int_t    GetMinITSHits()                     const {return fMinITSHits;}
352   //
353   Double_t GetUpdCalls()                       const {return fUpdCalls;}
354   void     InitMCWatchHistos();
355   TH2F*    GetHMCLrResidRPhi()                 const {return fHMCLrResidRPhi;}
356   TH2F*    GetHMCLrResidZ()                    const {return fHMCLrResidZ;}
357   TH2F*    GetHMCLrChi2()                      const {return fHMCLrChi2;}
358   //
359   void     PrintITS(Option_t* opt="") const {for (int i=0;i<=fLastActiveITSLayer;i++) if (!GetLayer(i)->IsDead()) GetLayer(i)->Print(opt);}
360   static void SetVtxConstraint(double d=-1, double z=-1) {fgVtxConstraint[0]=d; fgVtxConstraint[1]=z;}
361   //
362   void CalcHardSearchLimits(double dzv);
363   void SetMaxSeedToPropagate(Int_t n=300) {fMaxSeedToPropagate = n;}
364
365  protected:
366  
367   Int_t fNLayers;        // total number of layers in the model
368   Int_t fNActiveLayers;  // number of active layers in the model
369   Int_t fNActiveITSLayers;  // number of active ITS layers in the model
370   Int_t fLastActiveLayer;       // id of last active layer
371   Int_t fLastActiveITSLayer;    // id of last active ITS layer
372   Int_t fLastActiveLayerTracked;    // id of last active layer really used for tracking of given pt
373   TList fLayers;                // List of layer pointers
374   Float_t fBFieldG;             // Magnetic Field in Gauss (set in Tesla)
375   Float_t fLhcUPCscale;         // UltraPeripheralElectrons: scale from RHIC to LHC
376   Float_t fIntegrationTime;     // electronics integration time
377   Float_t fConfLevel;           // Confidence Level for the tracking
378   Float_t fAvgRapidity;         // rapidity of the track (= mean)
379   Float_t fParticleMass;        // Particle used for tracking. Standard: mass of pion
380   Double_t fMaxRadiusSlowDet;   // Maximum radius for slow detectors.  Fast detectors 
381                                 // and only fast detectors reside outside this radius.
382   Int_t fAtLeastCorr;     // min. number of correct hits for the track to be "good"
383   Int_t fAtLeastFake;     // min. number of fake hits for the track to be "fake"
384
385   Int_t fdNdEtaCent;       // Multiplicity
386   //
387   // reconstruction settings
388   Double_t fMaxChi2Cl;   // max cluster-track chi2 
389   Double_t fMaxNormChi2NDF;// max chi2/NDF to accept
390   Int_t    fMinITSHits;  // min ITS hits in track to accept
391   //
392   Double_t fMaxChi2ClSQ; // precalculated sqrt(chi2);
393   //
394   Int_t    fMaxSeedToPropagate; // take 1st fMaxSeedToPropagate seeds from each layer
395   // background settings
396   Bool_t   fUseBackground; // do we want to simulate background?
397   Double_t fBgYMin;      // min Y for current layer bg generation
398   Double_t fBgYMax;      // max Y for current layer bg generation
399   Double_t fBgZMin;      // min Z for current layer bg generation
400   Double_t fBgZMax;      // max Z for current layer bg generation
401   TArrayD  fBgYMinTr;    // min Y for each seed at current layer
402   TArrayD  fBgYMaxTr;    // max Y for each seed at current layer
403   TArrayD  fBgZMinTr;    // min Z for each seed at current layer
404   TArrayD  fBgZMaxTr;    // max Z for each seed at current layer
405   Int_t    fNBgLimits;   // depth of limits array
406   //
407   enum {kMaxNDetectors = 200};
408  
409   Double_t fTransMomenta[40];                          // array of transverse momenta
410   Double_t fMomentumRes[40];                           // array of momentum resolution
411   Double_t fResolutionRPhi[40];                        // array of rphi resolution
412   Double_t fResolutionZ[40];                           // array of z resolution
413   Double_t fDetPointRes[kMaxNDetectors][40];           // array of rphi resolution per layer
414   Double_t fDetPointZRes[kMaxNDetectors][40];          // array of z resolution per layer
415   Double_t fEfficiency[3][40];                         // efficiency for different particles
416   Double_t fFake[3][40];                               // fake prob for different particles
417   //
418   KMCProbe fProbe;
419   Double_t fDensFactorEta;                             // density scaling for non-0 eta (precalculated
420   //
421   Double_t fUpdCalls;                  // number of kalman updates
422   TH2F*  fHMCLrResidRPhi;              // Residials on lr, rphi
423   TH2F*  fHMCLrResidZ;                 // Residials on lr, rphi
424   TH2F*  fHMCLrChi2;                   // track to clusyer chi2 on each lr, rphi
425   TArrayI fPattITS;                     // bit pattern of each group of active layers where hit is required
426   //
427   static Double_t fgVtxConstraint[2];  // if both positive, the vertex is used as constraint (accounted in chi2 but not in update)
428   ClassDef(KMCDetector,1);
429 };
430
431 //____________________________________________________________________________
432 inline Bool_t KMCProbe::PropagateToCluster(KMCCluster* cl, double b)
433 {
434   // propagate track to cluster frame
435   if (!Rotate(cl->GetPhi()) || !PropagateTo(cl->GetX(),b)) {
436     AliDebug(2,Form("Failed to propager track to cluster at phi=%.3f X=%.3f",cl->GetPhi(),cl->GetX()));
437     if (AliLog::GetGlobalDebugLevel()>1) Print();
438     return kFALSE;
439   }
440   return kTRUE;
441 }
442
443 //////////////////////////////////////////////////////////////////////////////////////////////
444 //--------------------------------------------------------------------------------------------
445 //
446 // Class to collect summary of certrain track performance
447 // The track can be selected according to min and max ITS clusters to accept (total, fake and or correct)
448 // and according to presence of (any and/or correct)  or absence of fakes clusters in certain groups
449 // of ITS layers.
450
451 class KMCTrackSummary : public TNamed
452 {
453  public:
454   KMCTrackSummary(const char* name=0,const char* title=0,int nlr=0);
455   virtual ~KMCTrackSummary();
456   //
457   void Add(const KMCTrackSummary* src, double scl=1);
458
459   virtual void  Print(Option_t* option = "") const;
460   KMCTrackSummary* MakeCopy(const char* pref) const;
461   void   SetNamePrefix(const char* pref="");
462
463   Bool_t CheckTrack(KMCProbe* trc);
464   void   AddTrack(KMCProbe* trc);
465   //
466   TH1*   GetHMCChi2()               const {return fHMCChi2;}
467   TH1*   GetHMCSigDCARPhi()         const {return fHMCSigDCARPhi;}
468   TH1*   GetHMCSigDCAZ()            const {return fHMCSigDCAZ;}
469   TH1*   GetHMCSigPt()              const {return fHMCSigPt;}
470   TH2*   GetHMCNCl()                const {return fHMCNCl;}
471   TH1*   GetHMCFakePatt()           const {return fHMCFakePatt;}
472   //
473   Double_t GetCountAcc()            const {return fCountAcc;}
474   Double_t GetCountTot()            const {return fCountTot;}
475   Double_t GetEff()                 const {return fCountTot>0 ? fCountAcc/fCountTot : 0.;}
476   Double_t GetEffErr()              const;
477   KMCProbe& GetAnProbe()      const {return (KMCProbe&)fAnProbe;}
478   KMCProbe& GetRefProbe()     const {return (KMCProbe&)fRefProbe;}
479   void SetAnProbe(KMCProbe* pr)     {if (pr) fAnProbe = *pr;}
480   void SetRefProbe(KMCProbe* pr)    {if (pr) fRefProbe = *pr;}
481   //
482   void SetMinMaxClITS(Int_t mn=0,Int_t mx=999)     {fMinMaxClITS[0]=mn; fMinMaxClITS[1]=mx;}
483   void SetMinMaxClITSFake(Int_t mn=0,Int_t mx=999) {fMinMaxClITSFake[0]=mn; fMinMaxClITSFake[1]=mx;}
484   void SetMinMaxClITSCorr(Int_t mn=0,Int_t mx=999) {fMinMaxClITSCorr[0]=mn; fMinMaxClITSCorr[1]=mx;}
485   //
486   Double_t GetUpdCalls()           const {return fCountTot>0 ? fUpdCalls/fCountTot : 0;}
487   void     AddUpdCalls(double v)         {fUpdCalls += v;}
488   //
489   void AddPatternITS(Int_t patt)           {AddPattern(patt,fPattITS);}
490   void AddPatternITSFakeExcl(Int_t patt)   {AddPattern(patt,fPattITSFakeExcl);}
491   void AddPatternITSCorr(Int_t patt)       {AddPattern(patt,fPattITSCorr);}
492   static UInt_t Bits(Bool_t l0=0, Bool_t l1=0, Bool_t l2=0, Bool_t l3=0, Bool_t l4=0, Bool_t l5=0, Bool_t l6=0, Bool_t l7=0,Bool_t  l8=0, Bool_t l9=0,
493                      Bool_t l10=0,Bool_t l11=0,Bool_t l12=0,Bool_t l13=0,Bool_t l14=0,Bool_t l15=0,Bool_t l16=0,Bool_t l17=0,Bool_t l18=0,Bool_t l19=0,
494                      Bool_t l20=0,Bool_t l21=0,Bool_t l22=0,Bool_t l23=0,Bool_t l24=0,Bool_t l25=0,Bool_t l26=0,Bool_t l27=0,Bool_t l28=0,Bool_t l29=0,
495                      Bool_t l30=0,Bool_t l31=0);
496   //
497   // protected:
498   void   AddPattern(Int_t patt,TArrayI& arr) {int n=arr.GetSize(); arr.Set(n+1); arr[n]=patt;}
499   Bool_t CheckPattern(UInt_t patt, Int_t cond) {return ((UInt_t)cond)&patt;}
500   Int_t  OutOfRange(Int_t n, Int_t range[2])  {if (n<range[0]) return -1; if (n>range[1]) return 1; return 0;}
501   static void   PrependTNamed(TNamed* obj, const char* nm=0, const char* tit=0);
502   //
503  private:
504  KMCTrackSummary(const KMCTrackSummary& src) :TNamed(src) {}
505  protected:
506   //
507   Int_t   fNITSLayers;                  // number of its layers
508   // track selection conditions
509   Int_t   fMinMaxClITS[2];              // min and max N ITS clusters (total)
510   Int_t   fMinMaxClITSFake[2];          // min and max N ITS fake clusters  
511   Int_t   fMinMaxClITSCorr[2];          // min and max N ITS correct clusters
512   //
513   TArrayI fPattITS;                     // bit pattern of each group of layers where hit is required
514   TArrayI fPattITSFakeExcl;             // bit pattern of each group of layers where fake hit is not tolerated
515   TArrayI fPattITSCorr;                 // bit pattern of each group of layers where correc hit is required
516   //
517   TH1F*  fHMCChi2;                     // track chi2
518   TH1F*  fHMCSigDCARPhi;               // DCA distribution in RPhi from MC
519   TH1F*  fHMCSigDCAZ;                  // DCA distribution in Z for from MC
520   TH1F*  fHMCSigPt;                    // Pt difference distribution from MC
521   TH2F*  fHMCNCl;                      // number of clusters
522   TH1F*  fHMCFakePatt;                 // Fakes pattern
523   //
524   Double_t fCountAcc;                  // track counter of accepted tracks
525   Double_t fCountTot;                  // track counter of all tracks
526   Double_t fUpdCalls;                  // total update calls
527   //
528   KMCProbe fRefProbe;                  // reference probe
529   KMCProbe fAnProbe;                   // analital probe
530   static Int_t fgSumCounter;           // global counter of bins
531
532   ClassDef(KMCTrackSummary,1)
533 };
534
535 inline Double_t KMCTrackSummary::GetEffErr() const {
536   //
537   if (fCountTot<1) return 0;
538   double r = fCountAcc/fCountTot;
539   double err = r*(1.-r)/fCountTot;
540   return err>0 ? TMath::Sqrt(err) : 0;
541 }
542
543 inline KMCTrackSummary* KMCTrackSummary::MakeCopy(const char* pref) const {
544   // create a copy, prepending all names with prefix
545   KMCTrackSummary* sm = (KMCTrackSummary*) Clone(this->GetName()); 
546   sm->SetNamePrefix(pref); 
547   return sm;
548
549
550 #endif