AliPhysics  9c7580a (9c7580a)
AliTrackComparisonESD.cxx
Go to the documentation of this file.
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 
18 // //
19 // ANALYSIS task to perrorm TPC calibration //
20 
21 // //
23 #include "AliTrackComparisonESD.h"
24 #include "AliTrackComparison.h"
25 #include "TChain.h"
26 #include "AliESDEvent.h"
27 #include "AliESDfriend.h"
28 #include "AliESDVertex.h"
29 #include "AliESDtrack.h"
30 #include "AliESDfriendTrack.h"
31 #include "AliExternalTrackParam.h"
32 #include "AliTrackPointArray.h"
33 #include "AliESDtrackCuts.h"
34 #include "AliTracker.h"
35 #include "AliESDCaloCluster.h"
36 #include "AliESDInputHandler.h"
37 #include "AliAnalysisManager.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliCalorimeterUtils.h"
40 #include "AliESDCaloCells.h"
41 #include "TFile.h"
42 #include "TSystem.h"
43 #include "TTimeStamp.h"
44 #include "AliHMPIDParam.h"
45 //#include <TGeoHMatrix>
46 #include "AliGeomManager.h"
47 #include "AliGRPManager.h"
48 #include <iostream>
49 
50 using std::cout;
51 using std::endl;
52 
53 ClassImp(AliTrackComparisonESD)
54 
55 //________________________________________________________________________
57  :AliAnalysisTask(),
58  fESD(0),
59  fESDCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
60  fESDfriend(0),
61  fCurrentRun(-1),
62  fDebugOutputPath(""),
63  fOutput(0),
64  fEMCAL(0),
65  fHMPID(0),
66  fTOF(0),
67  fGeom(0),
68  fCutR(20),
69  fCaloUtil(0)
70 {
71  //
72  // default constructor
73  //
74  for(Int_t i=0; i<4; i++) fTransMatrix[i]=0;
75 
76 }
77 
78 //________________________________________________________________________
80  :AliAnalysisTask(name,""),
81  fESD(0),
82  fESDCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
83  fESDfriend(0),
84  fCurrentRun(-1),
85  fDebugOutputPath(""),
86  fOutput(0),
87  fEMCAL(0),
88  fHMPID(0),
89  fTOF(0),
90  fGeom(0),
91  fCutR(20),
92  fCaloUtil(0)
93 {
94  //
95  // Constructor
96  //
97  DefineInput(0, TChain::Class());
98  DefineOutput(0, TObjArray::Class());
99 
100  for(Int_t i=0; i<4; i++) fTransMatrix[i]=0;
101 }
102 
103 //________________________________________________________________________
105  //
106  // destructor
107  //
108  printf("AliTrackComparisonESD::~AliTrackComparisonESD");
109 
110  if(fEMCAL) delete fEMCAL; fEMCAL=0;
111  if(fHMPID) delete fHMPID; fHMPID=0;
112  if(fTOF) delete fTOF; fTOF=0;
113  if(fOutput) fOutput->Delete();
114  //if(fOutput) delete fOutput;
115  if(fCaloUtil) delete fCaloUtil; fCaloUtil=0;
116  for(Int_t i=0; i<4; i++)
117  {
118  if(fTransMatrix[i]) {delete fTransMatrix[i];fTransMatrix[i]=0;}
119  }
120 }
121 
122 //________________________________________________________________________
124  //
125  //
126  //
127  TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
128  if (!tree) {
129  //Printf("ERROR: Could not read chain from input slot 0");
130  }
131  else {
132  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
133  if (!esdH) {
134  //Printf("ERROR: Could not get ESDInputHandler");
135  }
136  else {
137  //esdH->SetReadFriends(kTRUE);
138  //esdH->SetActiveBranches("ESDfriend");
139  fESD = (AliESDEvent*)esdH->GetEvent();
140  //Printf("*** CONNECTED NEW EVENT ****");
141  }
142  }
143 
144 }
145 
146 //________________________________________________________________________
148  //
149  //
150  //
151  //OpenFile(0, "RECREATE");
152  TFile *ftmp = OpenFile(0);
153  if(!ftmp)AliError(Form("File %s not found!",ftmp->GetName()));
154  fOutput=new TObjArray(0);
155  fOutput->SetOwner(kTRUE);
156 
157  fEMCAL = new AliTrackComparison("EMCAL","EMCAL");
158  fEMCAL->SetLayerID(20);
159  fEMCAL->SetFillAll(kFALSE);
160  fEMCAL->Init();
161  fHMPID = new AliTrackComparison("HMPID","HMPID");
162  fHMPID->SetRangeDY(-5,5);
163  fHMPID->SetRangeDZ(-5,5);
164  fHMPID->SetLayerID(18);
165  fHMPID->SetFillAll(kFALSE);
166  fHMPID->Init();
167  fTOF = new AliTrackComparison("TOF","TOF");
168  fTOF->SetRangeDY(-5,5);
169  fTOF->SetRange1Pt(-1,1);
170  fTOF->SetLayerID(15);
171  fTOF->SetFillAll(kFALSE);
172  fTOF->Init();
173 
174  fOutput->Add(fEMCAL);
175  fOutput->Add(fHMPID);
176  fOutput->Add(fTOF);
177 
178  Double_t rotationMatrix[4][9] = {{-0.014587, -0.999892, -0.002031, 0.999892, -0.014591, 0.001979, -0.002009, -0.002002, 0.999996},
179  {-0.014587, 0.999892, 0.002031, 0.999892, 0.014591, -0.001979, -0.002009, 0.002002, -0.999996},
180  {-0.345864, -0.938278, -0.003412, 0.938276, -0.345874, 0.003010, -0.004004, -0.002161, 0.999990},
181  {-0.345864, 0.938278, 0.003412, 0.938276, 0.345874, -0.003010, -0.004004, 0.002161, -0.999990}};
182 
183  Double_t translationMatrix[4][3] = {{0.351659, 447.576446, 176.269742},
184  {1.062577, 446.893974, -173.728870},
185  {-154.213287, 419.306156, 176.753692},
186  {-153.018950, 418.623681, -173.243605}};
187  for(Int_t imx=0; imx<4; imx++)
188  {
189  fTransMatrix[imx] = new TGeoHMatrix();
190  fTransMatrix[imx]->SetRotation(rotationMatrix[imx]);
191  fTransMatrix[imx]->SetTranslation(translationMatrix[imx]);
192  fTransMatrix[imx]->Print();
193  }
194 
195 
197  InitCaloUtil();
198 
199 
200  PostData(0,fOutput);
201 }
202 
203 //________________________________________________________________________
205  //
206  // Setup Event
207  //
208  // check if something to be done
209 
210  if(!fESD)
211  return kFALSE;
212 
213  if (fCurrentRun == fESD->GetRunNumber())
214  return kTRUE;
215  else
216  fCurrentRun = fESD->GetRunNumber();
217 
218  return kTRUE;
219 }
220 
221 //________________________________________________________________________
223  //
224  // Exec function
225  // Loop over tracks and call Process function
226  if (!fESD) {
227  //Printf("ERROR: fESD not available");
228  return;
229  }
230 
231  if(!SetupEvent()) return;
232 
233 
234  fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
235  if (!fESDfriend) {
236  //Printf("ERROR: fESDfriend not available");
237  return;
238  }
239 
240 
241  if ( fESDfriend->GetNumberOfTracks() <=0 ) {
242  Printf("ERROR: fESDfriend Tracks not available");
243  return;
244  }
245 
246 
247  //Get primary vertex
248  const AliESDVertex *vertex = fESD->GetPrimaryVertex();
249  if(!vertex) AliError("No primary vertex found!\n");
250  Double_t vPos[3];
251  vertex->GetXYZ(vPos);
252  if(TMath::Abs(vPos[2])>7) return;
253 
254 
255  //Get EMCAL clusters and cells
256  TRefArray *clusters = new TRefArray();
257  Int_t nclusters = fESD->GetEMCALClusters(clusters);
258  AliESDCaloCells *cells = fESD->GetEMCALCells();
259  // RecalClusterPos(clusters,cells);
260 // Float_t pos[3];
261 // for(Int_t icl=0; icl<nclusters; icl++)
262 // {
263 // AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(icl);
264 // cluster->GetPosition(pos);
265 // printf("cluster %d pass00 pos: (%5.3f,%5.3f,%5.3f,%5.3f)\n",icl,pos[0],pos[1],pos[2],TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]));
266 // }
267 
268 
269 
270  //Loop over tracks
271  Int_t ntracks = fESD->GetNumberOfTracks();
272  for(Int_t itr=0; itr<ntracks; itr++)
273  {
274  AliESDtrack *track = fESD->GetTrack(itr);
275  if(!track || !fESDCuts->AcceptTrack(track)) continue;
276 
277  AliESDfriendTrack *friendTrack = (AliESDfriendTrack*)track->GetFriendTrack();
278  if(!friendTrack) continue;
279  //printf(" --- %d < %d || %p | %p -- %p \n", itr, fESDfriend->GetNumberOfTracks(), track, fESDfriend, friendTrack);
280  ProcessTOF(track,friendTrack,vPos);
281  if(nclusters>0)ProcessEMCAL(track,friendTrack,clusters,cells,vPos);
282  ProcessHMPID(track,friendTrack,vPos);
283  }//End of track loop
284 
285  delete clusters;
286  PostData(0,fOutput);
287 }
288 
289 //________________________________________________________________________
290 void AliTrackComparisonESD::ProcessTOF(AliESDtrack *track, AliESDfriendTrack *friendTrack, Double_t *vPos){
291  //
292  // Process TPC-TOF extrapolation
293  //
294 
295  //printf("Enter function!\n");
296  if (track->GetTOFsignal()<=0) return;
297  if (!friendTrack) return;
298  if (!friendTrack->GetTPCOut()) return;
299 
300  AliExternalTrackParam *pTPC = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
301  if(!pTPC) return;
302 
303  const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
304  if (!points) return;
305  Int_t npoints = points->GetNPoints();
306  if(npoints>1000) return; //the default value is more than 30000, why not -1???
307  AliTrackPoint point;
308  //
309  Int_t counter=0;
310  for (Int_t ipoint=0;ipoint<npoints;ipoint++){
311  if(!points->GetPoint(point,ipoint)) continue;
312  if(AliGeomManager::VolUIDToLayer(point.GetVolumeID())==AliGeomManager::kTOF)
313  {
314  counter++;
315  fTOF->AddTracks(pTPC,&point,track->GetMass(),track->P(),vPos);
316  }
317  }
318  //Printf("# of track points in TOF: %d!\n",counter);
319  return;
320 }
321 
322 //________________________________________________________________________
323 void AliTrackComparisonESD::ProcessEMCAL(AliESDtrack *track, AliESDfriendTrack *friendTrack, TRefArray *clusters, AliESDCaloCells *cells, Double_t *vPos){
324  if(clusters->GetEntriesFast()==0) return;
325 
326  Double_t rEMCal = 438;
327  Double_t tmp=fCutR;
328  Int_t clsIndex=-1;
329  TVector3 clsV,trkV;
330 
331  AliExternalTrackParam *pTPC = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
332  if(!pTPC) return;
333 
334  Double_t trPos[3];
335  Float_t clPos[3];
336  AliExternalTrackParam *pTest = new AliExternalTrackParam(*pTPC);
337  if(!AliTracker::PropagateTrackToBxByBz(pTest, rEMCal , track->GetMass(), 1 , kFALSE,0.99,-1)) return;
338  if(!pTest->GetXYZ(trPos)) return;
339 
340  Double_t trPhi = TMath::ATan2(trPos[1],trPos[0])*TMath::RadToDeg();
341  //printf("trPhi = %5.3f | eta = %5.3f\n",trPhi,pTest->Eta());
342  if(trPhi<80 || trPhi>120) return;
343  if(TMath::Abs(pTest->Eta())>0.7) return;
344 
345  AliExternalTrackParam *p0=0;
346  AliTrackPoint *p1=new AliTrackPoint();
347  Int_t nclusters = clusters->GetEntries();
348  for(Int_t icl=0; icl<nclusters; icl++)
349  {
350  AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(icl);
351  if(!cluster) continue;
352  if(fCaloUtil->ClusterContainsBadChannel("EMCAL",cluster->GetCellsAbsId(),cluster->GetNCells()) ) continue;
353  if(!fCaloUtil->CheckCellFiducialRegion(cluster,cells,NULL,0) ) continue;
354 
355  cluster->GetPosition(clPos);
356  p0 = pTPC;
357  clsV.SetXYZ(clPos[0],clPos[1],clPos[2]);
358  Double_t alpha = ((int)(clsV.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
359  clsV.RotateZ(-alpha);
360  p0->Rotate(alpha);
361  if(!AliTrackerBase::PropagateTrackToBxByBz(p0,clsV.X(), track->GetMass(), 1,kFALSE,0.99,-1)) continue;
362  trkV.SetXYZ(p0->GetX(),p0->GetY(),p0->GetZ());
363  Double_t dist = TMath::Sqrt( TMath::Power(clsV.X()-trkV.X(),2)+TMath::Power(clsV.Y()-trkV.Y(),2)+TMath::Power(clsV.Z()-trkV.Z(),2) );
364  if(dist<tmp)
365  {
366  tmp=dist;
367  clsIndex=icl;
368  }
369  }
370 
371  if(clsIndex==-1) return;
372  AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(clsIndex);
373  cluster->GetPosition(clPos);
374  //printf("cluster pos: (%5.3f,%5.3f,%5.3f,%5.3f)\n",clPos[0],clPos[1],clPos[2],TMath::Sqrt(clPos[0]*clPos[0]+clPos[1]*clPos[1]));
375  p1->SetXYZ(clPos[0],clPos[1],clPos[2],0);
376  //printf("Found EMCAL point!\n");
377  fEMCAL->AddTracks(pTPC,p1,track->GetMass(),track->P(),vPos);
378 
379 
380  delete pTest;
381  delete p1;
382  return;
383 }
384 
385 //________________________________________________________________________
386 void AliTrackComparisonESD::ProcessHMPID(AliESDtrack *track, AliESDfriendTrack *friendTrack, Double_t *vPos){
387  //
388  // Process TPC-TOF extrapolation
389  //
390  if (track->GetHMPIDsignal()<=0) return;
391 
392  AliExternalTrackParam *pTPC = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
393  if(!pTPC) return;
394 
395  Int_t q, nph, ch;
396  Float_t x, y;
397  track->GetHMPIDmip(x,y,q,nph);
398  Double_t pHmp[3]={0}, pHmp3=0;
399  if (track->GetOuterHmpPxPyPz(pHmp))
400  pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
401 
402  ch = track->GetHMPIDcluIdx()/1000000;
403 
404  AliHMPIDParam *pParam = AliHMPIDParam::Instance();
405  TVector3 vG = pParam->Lors2Mars(ch,x,y);
406 
407  AliTrackPoint *p1 = new AliTrackPoint();
408  p1->SetXYZ(vG.X(),vG.Y(),vG.Z());
409  //printf("Found HMPID point!\n");
410  fHMPID->AddTracks(pTPC,p1,track->GetMass(),pHmp3,vPos);
411  delete p1;
412  return;
413 }
414 
415 
416 //________________________________________________________________________
417 void AliTrackComparisonESD::RecalClusterPos(TRefArray *clusters, AliESDCaloCells *cells){
418  //
419  //
420  //
421  Double_t iLocal[3], iGlobal[3];
422  Float_t cPos[3];
423  //Float_t pos[3];
424  Int_t nclusters = clusters->GetEntries();
425  for(Int_t icl=0; icl<nclusters; icl++)
426  {
427  AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(icl);
428  UShort_t *absId = cluster->GetCellsAbsId();
429  Int_t nCells = cluster->GetNCells();
430  for(Int_t i=0;i<3;i++)cPos[i]=0;
431  Double_t wTot=0;
432  for(Int_t iCell=0; iCell<nCells; iCell++)
433  {
434  Double_t cellEnergy = cells->GetCellAmplitude(absId[iCell]);
435  Double_t dist = 1.31*(log(cluster->E())+4.82+0.5);
436  fGeom->RelPosCellInSModule(absId[iCell],dist,iLocal[0],iLocal[1],iLocal[2]);
437  //fGeom->GetGlobal(iLocal,iGlobal,fGeom->GetSuperModuleNumber(absId[iCell]));
438  Int_t sm = fGeom->GetSuperModuleNumber(absId[iCell]);
439  //matrix[sm]->Print();
440  //cout<<"sm = "<<sm<<endl;
441  fTransMatrix[sm]->LocalToMaster(iLocal,iGlobal);
442 
443  Double_t w = TMath::Max( 0., 4.5 + TMath::Log( cellEnergy / cluster->E() ));
444  if(w>0.0)
445  {
446  wTot += w;
447  for(Int_t i=0; i<3; i++ )
448  cPos[i] += (w*iGlobal[i]);
449  }
450  }//End of cell loop
451  if(wTot>0)
452  {
453  for(int i=0; i<3; i++ )
454  {
455  cPos[i] /= wTot;
456  cluster->SetPosition(cPos);
457  }
458  }
459  //cluster->GetPosition(pos);
460  //printf("cluster %d pass10 pos: (%5.3f,%5.3f,%5.3f)\n",icl,pos[0],pos[1],pos[2]);
461  }//End of cluster loop
462 }
463 
464 //________________________________________________________________________
466  //
467  // Terminate
468  //
469  AliInfo("AliTrackComparisonESD::Terminate()\n");
470 
471 }
472 
473 //________________________________________________________________________
475  //
476  // According description in AliAnalisysTask this method is call
477  // on the slaves before sending data
478  //
479  Terminate("slave");
480  if(!fDebugOutputPath.IsNull()) {
482  }
483 }
484 
485 //________________________________________________________________________
487  if(li) return 1;
488  else return 1;
489 }
490 
491 //________________________________________________________________________
493  //
494  // Analyze the content of the task
495  //
496 
497 }
498 
499 //________________________________________________________________________
501  //
502  //
503  //
504 }
505 
506 //________________________________________________________________________
508  cout<<"Initialize bad channel map!"<<endl;
511  //fCaloUtil->SetNumberOfCellsFromPHOSBorder(2);
514  // SM0
520  // SM1
535  // SM2
549  //SM3
551 
552  fCaloUtil->Print("");
553  cout<<"Done initialization!"<<endl;
554 }
double Double_t
Definition: External.C:58
AliTrackComparison * fTOF
Bool_t ClusterContainsBadChannel(Int_t calo, UShort_t *cellList, Int_t nCells)
long long Long64_t
Definition: External.C:43
TGeoHMatrix * fTransMatrix[4]
void SwitchOnNoFiducialBorderInEMCALEta0()
void SetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow, Double_t c=1)
Int_t fCurrentRun
current esd friend
virtual void Exec(Option_t *option)
AliESDtrackCuts * fESDCuts
current esd
virtual void ConnectInputData(Option_t *option)
void SetNumberOfCellsFromEMCALBorder(Int_t n)
void SetRange1Pt(Double_t lowBin1Pt, Double_t upBin1Pt)
void ProcessHMPID(AliESDtrack *track, AliESDfriendTrack *friendTrack, Double_t *vPos)
void ProcessEMCAL(AliESDtrack *track, AliESDfriendTrack *friendTrack, TRefArray *clusters, AliESDCaloCells *cells, Double_t *vPos)
virtual void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Int_t AddTracks(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass)
void SetLayerID(Double_t layerID)
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
void ProcessTOF(AliESDtrack *track, AliESDfriendTrack *friendTrack, Double_t *vPos)
void SetRangeDY(Double_t lowBinDY, Double_t upBinDY)
void RecalClusterPos(TRefArray *clusters, AliESDCaloCells *cells)
AliEMCALGeometry * fGeom
void SetRangeDZ(Double_t lowBinDZ, Double_t upBinDZ)
AliESDfriend * fESDfriend
esd track cuts
AliCalorimeterUtils * fCaloUtil
AliTrackComparison * fEMCAL
AliTrackComparison * fHMPID
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Class with utils specific to calorimeter clusters/cells.
virtual void Terminate(Option_t *option)
virtual Long64_t Merge(TCollection *li)
void SetFillAll(Bool_t fillAll)
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Bool_t CheckCellFiducialRegion(AliVCluster *cluster, AliVCaloCells *cells) const