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 "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
  48 ClassImp(AliTrackComparisonESD)
  49
  50 //________________________________________________________________________
  51 AliTrackComparisonESD::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 //________________________________________________________________________
  76 AliTrackComparisonESD::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 //________________________________________________________________________
 103 AliTrackComparisonESD::~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 //________________________________________________________________________
 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   PostData(0,fOutput);
 193 }
 194
 195 //________________________________________________________________________
 196 Bool_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 //________________________________________________________________________
 246 void 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 //________________________________________________________________________
 314 void 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 //________________________________________________________________________
 357 void 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 //________________________________________________________________________
 393 void 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 //________________________________________________________________________
 423 void 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 //________________________________________________________________________
 468 void AliTrackComparisonESD::Terminate(Option_t */*option*/) {
 469   //
 470   // Terminate
 471   //
 472   AliInfo("AliTrackComparisonESD::Terminate()\n");
 473
 474 }
 475
 476 //________________________________________________________________________
 477 void 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 //________________________________________________________________________
 489 Long64_t AliTrackComparisonESD::Merge(TCollection *li) {
 490   if(li) return 1;
 491   else return 1;
 492 }
 493
 494 //________________________________________________________________________
 495 void AliTrackComparisonESD::Analyze() {
 496   //
 497   // Analyze the content of the task
 498   //
 499
 500 }
 501
 502 //________________________________________________________________________
 503 void AliTrackComparisonESD::RegisterDebugOutput(){
 504   //
 505   //
 506   //
 507 }