Added TRD script for automatic trending + change permissions on FMD script
[u/mrichter/AliRoot.git] / PWGPP / 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 using std::cout;
50 using std::endl;
51
52 ClassImp(AliTrackComparisonESD)
53
54 //________________________________________________________________________
55 AliTrackComparisonESD::AliTrackComparisonESD()
56   :AliAnalysisTask(),
57    fESD(0),
58    fESDCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
59    fESDfriend(0),
60    fCurrentRun(-1),
61    fDebugOutputPath(""),
62    fOutput(0),
63    fEMCAL(0),
64    fHMPID(0),
65    fTOF(0),
66    fGeom(0),
67    fCutR(20),
68    fCaloUtil(0)
69 {
70   //
71   // default constructor
72   // 
73   for(Int_t i=0; i<4; i++) fTransMatrix[i]=0;
74   
75 }
76
77 //________________________________________________________________________
78 AliTrackComparisonESD::AliTrackComparisonESD(const char *name) 
79   :AliAnalysisTask(name,""),
80    fESD(0),
81    fESDCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
82    fESDfriend(0),
83    fCurrentRun(-1),
84    fDebugOutputPath(""),
85    fOutput(0),
86    fEMCAL(0),
87    fHMPID(0),
88    fTOF(0),
89    fGeom(0),
90    fCutR(20),
91    fCaloUtil(0)
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 //________________________________________________________________________
103 AliTrackComparisonESD::~AliTrackComparisonESD() {
104   //
105   // destructor
106   //
107   printf("AliTrackComparisonESD::~AliTrackComparisonESD");
108
109   if(fEMCAL)  delete fEMCAL;  fEMCAL=0;
110   if(fHMPID)  delete fHMPID;  fHMPID=0;
111   if(fTOF)    delete fTOF;    fTOF=0;
112   if(fOutput) fOutput->Delete();
113   //if(fOutput) delete fOutput;
114   if(fCaloUtil) delete fCaloUtil; fCaloUtil=0;
115   for(Int_t i=0; i<4; i++)
116     {
117       if(fTransMatrix[i]) {delete fTransMatrix[i];fTransMatrix[i]=0;}
118     }
119 }
120
121 //________________________________________________________________________
122 void 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 //________________________________________________________________________
146 void 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
195   fCaloUtil = new AliCalorimeterUtils();
196   InitCaloUtil();
197
198
199   PostData(0,fOutput);
200 }
201
202 //________________________________________________________________________
203 Bool_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
217   return kTRUE;
218 }
219
220 //________________________________________________________________________
221 void 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 ) {
241     Printf("ERROR: fESDfriend Tracks not available");
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();
258   //  RecalClusterPos(clusters,cells);
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);
280       if(nclusters>0)ProcessEMCAL(track,friendTrack,clusters,cells,vPos);
281       ProcessHMPID(track,friendTrack,vPos);
282     }//End of track loop
283
284   delete clusters;
285   PostData(0,fOutput);
286 }
287
288 //________________________________________________________________________
289 void 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++){
310     if(!points->GetPoint(point,ipoint)) continue;
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);
318   return;
319 }
320
321 //________________________________________________________________________
322 void AliTrackComparisonESD::ProcessEMCAL(AliESDtrack *track, AliESDfriendTrack *friendTrack, TRefArray *clusters, AliESDCaloCells *cells, Double_t *vPos){
323   if(clusters->GetEntriesFast()==0) return;
324
325   Double_t rEMCal = 438;
326   Double_t tmp=fCutR;
327   Int_t clsIndex=-1;
328   TVector3 clsV,trkV;
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
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
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;
351       if(fCaloUtil->ClusterContainsBadChannel("EMCAL",cluster->GetCellsAbsId(),cluster->GetNCells()) ) continue;
352       if(!fCaloUtil->CheckCellFiducialRegion(cluster,cells,NULL,0) ) continue;
353
354       cluster->GetPosition(clPos);
355       p0 = pTPC;
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         }
368     }
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
378
379   delete pTest;
380   delete p1;
381   return;
382 }
383
384 //________________________________________________________________________
385 void 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;
411   return;
412 }
413
414
415 //________________________________________________________________________
416 void AliTrackComparisonESD::RecalClusterPos(TRefArray *clusters, AliESDCaloCells *cells){
417   //
418   //
419   //
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 //________________________________________________________________________
464 void AliTrackComparisonESD::Terminate(Option_t */*option*/) {
465   //
466   // Terminate
467   //
468   AliInfo("AliTrackComparisonESD::Terminate()\n");
469   
470 }
471
472 //________________________________________________________________________
473 void 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 //________________________________________________________________________
485 Long64_t AliTrackComparisonESD::Merge(TCollection *li) {
486   if(li) return 1;
487   else return 1;
488 }
489
490 //________________________________________________________________________
491 void AliTrackComparisonESD::Analyze() {
492   //
493   // Analyze the content of the task
494   //
495
496 }
497
498 //________________________________________________________________________
499 void AliTrackComparisonESD::RegisterDebugOutput(){
500   //
501   //
502   //
503 }
504
505 //________________________________________________________________________
506 void 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 }