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