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