Moving some functionality to selection base class.
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtSelectorPhos.cxx
1 //_________________________________________________________________________
2 //  Utility Class for transverse energy studies
3 //  Selector Base class for PHOS
4 //  -  
5 // implementation file
6 //
7 //*-- Authors: Oystein Djuvsland (Bergen)
8 //_________________________________________________________________________
9 #include "AliAnalysisEtSelectorPhos.h"
10 #include "AliAnalysisEtCuts.h"
11 #include "AliESDCaloCluster.h"
12 #include "AliESDEvent.h"
13 #include "TRefArray.h"
14 #include "AliPHOSGeometry.h"
15 #include "TH2I.h"
16 #include "TFile.h"
17 #include "TMath.h"
18 #include "TParticle.h"
19 #include "AliLog.h"
20 #include <iostream>
21
22 ClassImp(AliAnalysisEtSelectorPhos)
23
24 AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
25 ,fGeoUtils(0)
26 ,fBadMapM2(0)
27 ,fBadMapM3(0)
28 ,fBadMapM4(0)
29 ,fMatrixInitialized(kFALSE)
30 {
31   
32 }
33
34 AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
35 {
36
37 }
38
39 TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
40 { // Get clusters
41   if(!fClusterArray) fClusterArray = new TRefArray;
42   
43   if(fClusterArray)
44   {
45     fEvent->GetPHOSClusters(fClusterArray);
46   }
47   else
48   {
49     Printf("Could not initialize cluster array");
50   }
51   
52   return fClusterArray;
53 }
54
55 Int_t AliAnalysisEtSelectorPhos::Init(const AliESDEvent* event)
56 { // Init
57   
58   AliAnalysisEtSelector::Init(event);
59   Printf("Initializing selector for run: %d", event->GetRunNumber());
60   int res = LoadGeometry();
61   if(res) return -1;
62   if(LoadBadMaps()) return -1;
63   fInitialized = kTRUE;
64   if (!fMatrixInitialized)
65     {
66         Printf("INITIALIZING MISALIGNMENT MATRICES");
67         for (Int_t mod=0; mod<5; mod++) {
68             
69             if (!event->GetPHOSMatrix(mod))
70             {
71               Printf("Could not find geo matrix for module %d", mod);
72               continue;
73             }
74             fMatrixInitialized = kTRUE;
75             fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
76             Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
77         }
78     }
79   return 0;
80 }
81
82 Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const AliESDCaloCluster& cluster) const
83 {
84   
85 //    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
86   return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
87 }
88
89 Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const TParticle& part) const
90 {
91 //    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
92     return part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
93 }
94
95
96 Bool_t AliAnalysisEtSelectorPhos::PassDistanceToBadChannelCut(const AliESDCaloCluster& cluster) const
97 { // cut distance to bad channel
98   if(!fMatrixInitialized)
99   {
100     Printf("Misalignment matrices are not initialized");
101     return kFALSE;
102   }
103     Float_t gPos[3];
104     cluster.GetPosition(gPos);
105     Int_t relId[4];
106     TVector3 glVec(gPos);
107     fGeoUtils->GlobalPos2RelId(glVec, relId);
108
109     //std::cout << "In phos distance to bad channel cut!" << std::endl;
110     TVector3 locVec;
111     fGeoUtils->Global2Local(locVec, glVec, relId[0]);
112 //    std::cout << fGeoUtils << std::endl;
113     //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
114     //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
115     for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
116     {
117         for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
118         {
119             if (relId[0] == 3)
120             {
121                 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
122                 {
123                     Int_t tmpRel[4];
124                     tmpRel[0] = 3;
125                     tmpRel[1] = 0;
126                     tmpRel[2] = x+1;
127                     tmpRel[3] = z+1;
128                     
129                     Float_t tmpX;
130                     Float_t tmpZ;
131                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
132
133                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
134                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
135                     
136                     if (distance < fCuts->GetPhosBadDistanceCut())
137                     {
138 //                    std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
139                         return kFALSE;
140                     }
141                 }
142             }
143             if (relId[0] == 2)
144             {
145                 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
146                 {
147                     Int_t tmpRel[4];
148                     tmpRel[0] = 2;
149                     tmpRel[1] = 0;
150                     tmpRel[2] = x+1;
151                     tmpRel[3] = z+1;
152                     
153                     Float_t tmpX;
154                     Float_t tmpZ;
155                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
156
157                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
158
159 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
160                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
161                     if (distance < fCuts->GetPhosBadDistanceCut())
162                     {
163 //                    std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
164                         return kFALSE;
165                     }
166                 }
167             }
168             if (relId[0] == 1)
169             {
170                 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
171                 {
172                     Int_t tmpRel[4];
173                     tmpRel[0] = 1;
174                     tmpRel[1] = 0;
175                     tmpRel[2] = x+1;
176                     tmpRel[3] = z+1;
177                     
178                     Float_t tmpX;
179                     Float_t tmpZ;
180                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
181
182                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
183
184 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
185                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
186                     if (distance < fCuts->GetPhosBadDistanceCut())
187                     {
188 //                      std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
189                         return kFALSE;
190                     }
191                 }
192             }
193
194         }
195     }
196
197     return kTRUE;
198
199 }
200
201 Bool_t AliAnalysisEtSelectorPhos::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
202 { // cut track matching
203
204   if(!fMatrixInitialized)
205   {
206     Printf("Misalignment matrices are not initialized");
207     return kFALSE;
208   }
209   
210   // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
211   
212   Int_t nTracksMatched = cluster.GetNTracksMatched();
213   if(nTracksMatched == 0) return kTRUE;
214   
215   Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
216   if(trackMatchedIndex < 0) return kTRUE;
217   
218   AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
219   Double_t dx = cluster.GetTrackDx();
220   Double_t dz = cluster.GetTrackDz();
221   Double_t pt = track->Pt();
222   Int_t charge = track->Charge();
223   
224   Double_t meanX=0;
225   Double_t meanZ=0.;
226   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
227               6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
228   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
229   
230   Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
231
232   if(mf<0.){ //field --
233     meanZ = -0.468318 ;
234     if(charge>0)
235       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
236     else
237       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
238   }
239   else{ //Field ++
240     meanZ= -0.468318;
241     if(charge>0)
242       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
243     else
244       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;     
245   }
246
247   Double_t rz=(dz-meanZ)/sz ;
248   Double_t rx=(dx-meanX)/sx ;
249   Double_t r = TMath::Sqrt(rx*rx+rz*rz);
250   if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
251   
252   return kTRUE;
253  
254 }
255
256 int AliAnalysisEtSelectorPhos::LoadGeometry()
257 { // load geometry
258
259   fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
260     // ifstream f("badchannels.txt", ios::in);
261   return 0;
262 }
263
264 int AliAnalysisEtSelectorPhos::LoadBadMaps()
265 { // load bad maps
266 TFile *f = TFile::Open("badchannels.root", "READ");
267
268     if(!f)
269     {
270       std::cout << "Could not open badchannels.root" << std::endl;
271       return -1;
272     }
273
274     fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
275     if(!fBadMapM2) 
276     {
277       std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
278     }
279     fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
280     if(!fBadMapM2) 
281     {
282       std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
283     }
284     
285     fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
286     if(!fBadMapM4) 
287     {
288       std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
289     }
290     
291     
292     return 0;
293     
294 }
295
296
297 Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part) const
298 {
299   return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut() 
300           && part.Phi() < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
301           && part.Phi() > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
302 }
303
304 Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track) const
305 {
306   return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut() &&
307            track.Phi() > fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180. &&
308            track.Phi() < fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
309 }
310