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