o First Version of TRDnSigma implementation (Xianguo) o still requires some catching...
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtSelectorPhos.cxx
CommitLineData
4d376d01 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//_________________________________________________________________________
ef647350 9#include "AliAnalysisEtSelectorPhos.h"
10#include "AliAnalysisEtCuts.h"
11#include "AliESDCaloCluster.h"
f61cec2f 12#include "AliESDEvent.h"
ef647350 13#include "TRefArray.h"
14#include "AliPHOSGeometry.h"
15#include "TH2I.h"
16#include "TFile.h"
17#include "TMath.h"
f61cec2f 18#include "TParticle.h"
ef647350 19#include "AliLog.h"
20#include <iostream>
32503dac 21#include "AliStack.h"
f61cec2f 22
ef647350 23ClassImp(AliAnalysisEtSelectorPhos)
f61cec2f 24
ef647350 25AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
26,fGeoUtils(0)
27,fBadMapM2(0)
28,fBadMapM3(0)
29,fBadMapM4(0)
f61cec2f 30,fMatrixInitialized(kFALSE)
c31562f7 31{
32
33}
34
35AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(): AliAnalysisEtSelector()
36,fGeoUtils(0)
37,fBadMapM2(0)
38,fBadMapM3(0)
39,fBadMapM4(0)
40,fMatrixInitialized(kFALSE)
ef647350 41{
42
43}
44
45AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
46{
47
48}
49
50TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
4d376d01 51{ // Get clusters
f61cec2f 52 if(!fClusterArray) fClusterArray = new TRefArray;
53
54 if(fClusterArray)
55 {
56 fEvent->GetPHOSClusters(fClusterArray);
57 }
58 else
59 {
60 Printf("Could not initialize cluster array");
61 }
62
63 return fClusterArray;
ef647350 64}
65
f61cec2f 66Int_t AliAnalysisEtSelectorPhos::Init(const AliESDEvent* event)
4d376d01 67{ // Init
ef647350 68
f61cec2f 69 AliAnalysisEtSelector::Init(event);
70 Printf("Initializing selector for run: %d", event->GetRunNumber());
ef647350 71 int res = LoadGeometry();
72 if(res) return -1;
f61cec2f 73 if(LoadBadMaps()) return -1;
74 fInitialized = kTRUE;
75 if (!fMatrixInitialized)
76 {
77 Printf("INITIALIZING MISALIGNMENT MATRICES");
78 for (Int_t mod=0; mod<5; mod++) {
79
80 if (!event->GetPHOSMatrix(mod))
81 {
82 Printf("Could not find geo matrix for module %d", mod);
83 continue;
84 }
85 fMatrixInitialized = kTRUE;
86 fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
87 Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
88 }
89 }
90 return 0;
ef647350 91}
92
86e7d5db 93Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const AliESDCaloCluster& cluster) const
ef647350 94{
b2c10007 95
cfe23ff0 96 Float_t pos[3];
97 cluster.GetPosition(pos);
98 TVector3 cp(pos);
b2c10007 99// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
cfe23ff0 100 return TMath::Sin(cp.Theta())*cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
ef647350 101}
102
86e7d5db 103Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const TParticle& part) const
f61cec2f 104{
b2c10007 105// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
cfe23ff0 106 return TMath::Sin(part.Theta())*part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
f61cec2f 107}
108
02d47689 109Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(Double_t e) const
110{
111 return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
112}
113
f61cec2f 114
86e7d5db 115Bool_t AliAnalysisEtSelectorPhos::PassDistanceToBadChannelCut(const AliESDCaloCluster& cluster) const
4d376d01 116{ // cut distance to bad channel
f61cec2f 117 if(!fMatrixInitialized)
118 {
119 Printf("Misalignment matrices are not initialized");
120 return kFALSE;
121 }
ef647350 122 Float_t gPos[3];
123 cluster.GetPosition(gPos);
124 Int_t relId[4];
125 TVector3 glVec(gPos);
126 fGeoUtils->GlobalPos2RelId(glVec, relId);
127
b2c10007 128 //std::cout << "In phos distance to bad channel cut!" << std::endl;
ef647350 129 TVector3 locVec;
130 fGeoUtils->Global2Local(locVec, glVec, relId[0]);
131// std::cout << fGeoUtils << std::endl;
132 //std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
133 //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
134 for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
135 {
136 for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
137 {
138 if (relId[0] == 3)
139 {
140 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
141 {
142 Int_t tmpRel[4];
143 tmpRel[0] = 3;
144 tmpRel[1] = 0;
145 tmpRel[2] = x+1;
146 tmpRel[3] = z+1;
147
148 Float_t tmpX;
149 Float_t tmpZ;
150 fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
151
152 Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
153 //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
154
155 if (distance < fCuts->GetPhosBadDistanceCut())
156 {
157// std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
158 return kFALSE;
159 }
160 }
161 }
162 if (relId[0] == 2)
163 {
164 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
165 {
166 Int_t tmpRel[4];
167 tmpRel[0] = 2;
168 tmpRel[1] = 0;
169 tmpRel[2] = x+1;
170 tmpRel[3] = z+1;
171
172 Float_t tmpX;
173 Float_t tmpZ;
174 fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
175
176 Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
177
178// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
179 //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
180 if (distance < fCuts->GetPhosBadDistanceCut())
181 {
182// std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
183 return kFALSE;
184 }
185 }
186 }
187 if (relId[0] == 1)
188 {
189 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
190 {
191 Int_t tmpRel[4];
192 tmpRel[0] = 1;
193 tmpRel[1] = 0;
194 tmpRel[2] = x+1;
195 tmpRel[3] = z+1;
196
197 Float_t tmpX;
198 Float_t tmpZ;
199 fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
200
201 Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
202
203// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
204 //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
205 if (distance < fCuts->GetPhosBadDistanceCut())
206 {
207// std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
208 return kFALSE;
209 }
210 }
211 }
212
213 }
214 }
215
216 return kTRUE;
217
218}
219
86e7d5db 220Bool_t AliAnalysisEtSelectorPhos::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
4d376d01 221{ // cut track matching
ef647350 222
f61cec2f 223 if(!fMatrixInitialized)
224 {
225 Printf("Misalignment matrices are not initialized");
226 return kFALSE;
227 }
ef647350 228
229 // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
6a152780 230
ef647350 231 Int_t nTracksMatched = cluster.GetNTracksMatched();
b2c10007 232 if(nTracksMatched == 0) return kTRUE;
ef647350 233
234 Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
235 if(trackMatchedIndex < 0) return kTRUE;
236
237 AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
ea3ce170 238 if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random. It takes ~460 MeV to reach the EMCal
ef647350 239 Double_t dx = cluster.GetTrackDx();
240 Double_t dz = cluster.GetTrackDz();
241 Double_t pt = track->Pt();
242 Int_t charge = track->Charge();
243
244 Double_t meanX=0;
245 Double_t meanZ=0.;
246 Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
247 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);
248 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) ;
249
250 Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
251
252 if(mf<0.){ //field --
253 meanZ = -0.468318 ;
254 if(charge>0)
255 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)) ;
256 else
257 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)) ;
258 }
259 else{ //Field ++
260 meanZ= -0.468318;
261 if(charge>0)
262 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)) ;
263 else
264 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)) ;
265 }
266
267 Double_t rz=(dz-meanZ)/sz ;
268 Double_t rx=(dx-meanX)/sx ;
f61cec2f 269 Double_t r = TMath::Sqrt(rx*rx+rz*rz);
ef647350 270 if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
271
272 return kTRUE;
273
274}
275
276int AliAnalysisEtSelectorPhos::LoadGeometry()
4d376d01 277{ // load geometry
ef647350 278
279 fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
280 // ifstream f("badchannels.txt", ios::in);
281 return 0;
282}
283
284int AliAnalysisEtSelectorPhos::LoadBadMaps()
4d376d01 285{ // load bad maps
ef647350 286TFile *f = TFile::Open("badchannels.root", "READ");
287
288 if(!f)
289 {
290 std::cout << "Could not open badchannels.root" << std::endl;
291 return -1;
292 }
f61cec2f 293
ef647350 294 fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
f61cec2f 295 if(!fBadMapM2)
ef647350 296 {
297 std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
298 }
299 fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
f61cec2f 300 if(!fBadMapM2)
ef647350 301 {
302 std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
303 }
304
305 fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
f61cec2f 306 if(!fBadMapM4)
ef647350 307 {
308 std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
309 }
ef647350 310
311
f61cec2f 312 return 0;
ef647350 313
314}
f61cec2f 315
f61cec2f 316
43dd5a38 317Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part)
f61cec2f 318{
43dd5a38 319 float myphi = part.Phi();
320 myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
f61cec2f 321 return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
43dd5a38 322 && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
323 && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
f61cec2f 324}
325
43dd5a38 326Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track)
f61cec2f 327{
43dd5a38 328 float myphi = track.Phi();
329 myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
330 return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
331 && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
332 && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
f61cec2f 333}
b2c10007 334
31c813d5 335
43dd5a38 336Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
31c813d5 337{
338 Float_t pos[3];
339 cluster.GetPosition(pos);
340 TVector3 cp(pos);
43dd5a38 341 float myphi = cp.Phi();
342 myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
343 return TMath::Abs(cp.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
344 && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
345 && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
31c813d5 346}
32503dac 347UInt_t AliAnalysisEtSelectorPhos::GetLabel(const AliESDCaloCluster *cluster, AliStack *stack){
348 //Finds primary and estimates if it unique one?
349 Int_t n=cluster->GetNLabels() ;
32503dac 350 Int_t iMax=0;
beb92504 351
352 if (n > 0) {
353 Double_t* Ekin= new Double_t[n] ;
354 for(Int_t i=0; i<n; i++){
355 TParticle* p= stack->Particle(cluster->GetLabelAt(i)) ;
356 Ekin[i]=p->P() ; // estimate of kinetic energy
357 if(p->GetPdgCode()==-2212 || p->GetPdgCode()==-2112){
358 Ekin[i]+=1.8 ; //due to annihilation
359 }
360 }
361 Double_t eMax=0.;//eSubMax=0. ;
362 for(Int_t i=0; i<n; i++){
363 if(Ekin[i]>eMax){
364 // eSubMax=eMax;
365 eMax=Ekin[i];
366 iMax=i;
367 }
32503dac 368 }
32503dac 369// if(eSubMax>0.8*eMax)//not obvious primary
370// sure=kFALSE;
371// else
372// sure=kTRUE;
18875e70 373 delete [] Ekin;
beb92504 374
375 } // n>0
6844c491 376 UInt_t correctLabel = cluster->GetLabelAt(iMax);
377 correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);
378// //Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
379// if( stack->IsSecondaryFromMaterial(correctLabel) && correctLabel>0){//if this is flagged as a secondary then we look to see where it really came from
380// TParticle *hitParticle = stack->Particle(correctLabel);
381// if(hitParticle){
382// Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >400);
383// if(partVtxSecondary){//at this point we have something which converted near the detector. Let's find the mother particle
384// UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
385// if(mothIdx>0){
386// TParticle *mother = stack->Particle(mothIdx);
387// if(mother){
388// partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >400);
389// if(!partVtxSecondary) return mothIdx;
390// else{
391// }
392 // if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
393 // if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
394 // //cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
395 // //cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
396 // //so now we know that the particle originated near the cover and within the acceptance of the detector
397 // return mothIdx;
398 // }
399 // }
400// }
401// }
402// }
403
404// }
405// }
406 return correctLabel ; // DS: should this line be inside n>0 check, and return another value if n<=0 ?
32503dac 407
408}