[u/mrichter/AliRoot.git] / PWG1 / AliTrackComparisonESD.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// ANALYSIS task to perrorm TPC calibration                                  //

//                                                                           //
///////////////////////////////////////////////////////////////////////////////
#include "AliTrackComparisonESD.h"
#include "AliTrackComparison.h"
#include "TChain.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliESDfriendTrack.h"
#include "AliExternalTrackParam.h"
#include "AliTrackPointArray.h"
#include "AliESDtrackCuts.h"
#include "AliTracker.h"
#include "AliESDCaloCluster.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisManager.h"
#include "AliEMCALGeometry.h"
#include "TFile.h"
#include "TSystem.h"
#include "TTimeStamp.h"
#include "AliHMPIDParam.h"
//#include <TGeoHMatrix>
#include "AliGeomManager.h"
//#include "AliCDBManager.h"
//#include "AliGRPManager.h"

ClassImp(AliTrackComparisonESD)

//________________________________________________________________________
AliTrackComparisonESD::AliTrackComparisonESD()
  :AliAnalysisTask(),
   fESD(0),
   fESDCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
   fESDfriend(0),
   fCurrentRun(-1),
   fDebugOutputPath(""),
   //   fOcdbPath("local:///lustre/alice/alien/alice/data/2010/OCDB"),
   fOutput(0),
   fEMCAL(0),
   fHMPID(0),
   fTOF(0),
   fGeom(0),
   fCutX(10),
   fCutY(10),
   fCutZ(10)
{
  //
  // default constructor
  // 
  for(Int_t i=0; i<4; i++) fTransMatrix[i]=0;
  
}

//________________________________________________________________________
AliTrackComparisonESD::AliTrackComparisonESD(const char *name) 
  :AliAnalysisTask(name,""),
   fESD(0),
   fESDCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
   fESDfriend(0),
   fCurrentRun(-1),
   fDebugOutputPath(""),
   //   fOcdbPath("local:///lustre/alice/alien/alice/data/2010/OCDB"),
   fOutput(0),
   fEMCAL(0),
   fHMPID(0),
   fTOF(0),
   fGeom(0),
   fCutX(10),
   fCutY(10),
   fCutZ(10)
{
  //
  // Constructor
  //
  DefineInput(0, TChain::Class());
  DefineOutput(0, TObjArray::Class());

  for(Int_t i=0; i<4; i++) fTransMatrix[i]=0;
}

//________________________________________________________________________
AliTrackComparisonESD::~AliTrackComparisonESD() {
  //
  // destructor
  //
  printf("AliTrackComparisonESD::~AliTrackComparisonESD");
  if(fOutput) delete fOutput; fOutput=0;
  if(fEMCAL)  delete fEMCAL;  fEMCAL=0;
  if(fHMPID)  delete fHMPID;  fHMPID=0;
  if(fTOF)    delete fTOF;    fTOF=0;
  for(Int_t i=0; i<4; i++)
    {
      if(fTransMatrix[i]) {delete fTransMatrix[i];fTransMatrix[i]=0;}
    }
}

//________________________________________________________________________
void AliTrackComparisonESD::ConnectInputData(Option_t *) {
  //
  //
  //
  TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
  if (!tree) {
    //Printf("ERROR: Could not read chain from input slot 0");
  } 
  else {
    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
    if (!esdH) {
      //Printf("ERROR: Could not get ESDInputHandler");
    } 
    else {
      //esdH->SetReadFriends(kTRUE);
      //esdH->SetActiveBranches("ESDfriend");
      fESD = esdH->GetEvent();
      //Printf("*** CONNECTED NEW EVENT ****");
    }
  }

}

//________________________________________________________________________
void AliTrackComparisonESD::CreateOutputObjects() {
  //
  //
  //
  //OpenFile(0, "RECREATE");
  TFile *ftmp = OpenFile(0);
  if(!ftmp)AliError(Form("File %s not found!",ftmp->GetName()));
  fOutput=new TObjArray(0);
  fOutput->SetOwner(kTRUE);

  fEMCAL = new AliTrackComparison("EMCAL","EMCAL");
  fEMCAL->SetLayerID(20);
  fEMCAL->SetFillAll(kFALSE);
  fEMCAL->Init();
  fHMPID = new AliTrackComparison("HMPID","HMPID");
  fHMPID->SetRangeDY(-5,5);
  fHMPID->SetRangeDZ(-5,5);
  fHMPID->SetLayerID(18);
  fHMPID->SetFillAll(kFALSE);
  fHMPID->Init();
  fTOF   = new AliTrackComparison("TOF","TOF");
  fTOF->SetRangeDY(-5,5);
  fTOF->SetRange1Pt(-1,1);
  fTOF->SetLayerID(15);
  fTOF->SetFillAll(kFALSE);
  fTOF->Init();

  fOutput->Add(fEMCAL);
  fOutput->Add(fHMPID);
  fOutput->Add(fTOF);

  Double_t rotationMatrix[4][9] = {{-0.014587, -0.999892, -0.002031, 0.999892, -0.014591,  0.001979, -0.002009, -0.002002,  0.999996},
				   {-0.014587,  0.999892,  0.002031, 0.999892,  0.014591, -0.001979, -0.002009,  0.002002, -0.999996},
				   {-0.345864, -0.938278, -0.003412, 0.938276, -0.345874,  0.003010, -0.004004, -0.002161,  0.999990},
				   {-0.345864,  0.938278,  0.003412, 0.938276,  0.345874, -0.003010, -0.004004,  0.002161, -0.999990}};

  Double_t translationMatrix[4][3] = {{0.351659, 447.576446,  176.269742},
				      {1.062577, 446.893974, -173.728870},
				      {-154.213287, 419.306156,  176.753692},
				      {-153.018950, 418.623681, -173.243605}};
  for(Int_t imx=0; imx<4; imx++)
    {
      fTransMatrix[imx] = new TGeoHMatrix();
      fTransMatrix[imx]->SetRotation(rotationMatrix[imx]);
      fTransMatrix[imx]->SetTranslation(translationMatrix[imx]);
      fTransMatrix[imx]->Print();
    }


  PostData(0,fOutput);
}

//________________________________________________________________________
Bool_t AliTrackComparisonESD::SetupEvent() {
  //
  // Setup Event
  //
  // check if something to be done

  if(!fESD)
    return kFALSE;

  if (fCurrentRun == fESD->GetRunNumber())
    return kTRUE;
  else
    fCurrentRun = fESD->GetRunNumber();

//   // OCDB
//   printf("setting run to %d\n",fCurrentRun);
//   AliCDBManager::Instance()->SetDefaultStorage(fOcdbPath.Data());
//   AliCDBManager::Instance()->SetRun(fCurrentRun); 

//   // magnetic field
//   if ( !TGeoGlobalMagField::Instance()->GetField() ) {
//     printf("Loading field map...\n");
//     AliGRPManager grpMan;
//     if( !grpMan.ReadGRPEntry() ) { 
//       printf("Cannot get GRP entry\n"); 
//       return kFALSE;
//     }
//     if( !grpMan.SetMagField() ) { 
//       printf("Problem with magnetic field setup\n"); 
//       return kFALSE;
//     }
//   }

//   // geometry
//   printf("Loading geometry...\n");
//   AliGeomManager::LoadGeometry();
//   if( !AliGeomManager::ApplyAlignObjsFromCDB("GRP ITS TPC TRD TOF PHOS EMCAL HMPID") ) {
//     //printf("Problem with align objects\n");
//   }
  fGeom =  AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
  printf("%s\n",fGeom->GetName());
  if(!fGeom)
    {
      printf("EMCAL geometry not loaded!\n");
      return kFALSE;
    }
  return kTRUE;
}

//________________________________________________________________________
void AliTrackComparisonESD::Exec(Option_t *) {
  //
  // Exec function
  // Loop over tracks and call  Process function
  if (!fESD) {
    //Printf("ERROR: fESD not available");
    return;
  }

  if(!SetupEvent()) return;


  fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
  if (!fESDfriend) {
    //Printf("ERROR: fESDfriend not available");
    return;
  }


  if ( fESDfriend->GetNumberOfTracks() <=0 ) {
    //Printf("ERROR: fESDfriend Tracks not available");
    return;
  }


  //Get primary vertex
  const AliESDVertex *vertex = fESD->GetPrimaryVertex();
  if(!vertex) AliError("No primary vertex found!\n");
  Double_t vPos[3];
  vertex->GetXYZ(vPos);
  if(TMath::Abs(vPos[2])>7) return;

  
  //Get EMCAL clusters and cells
  TRefArray *clusters = new TRefArray();
  Int_t nclusters = fESD->GetEMCALClusters(clusters);
  AliESDCaloCells *cells = fESD->GetEMCALCells();
  RecalClusterPos(clusters,cells);
//   Float_t pos[3];
//   for(Int_t icl=0; icl<nclusters; icl++)
//     {
//       AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(icl);
//       cluster->GetPosition(pos);
//       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]));
//     }


  //Loop over tracks
  Int_t ntracks = fESD->GetNumberOfTracks();
  for(Int_t itr=0; itr<ntracks; itr++)
    {
      AliESDtrack *track = fESD->GetTrack(itr);
      if(!track || !fESDCuts->AcceptTrack(track)) continue;

      AliESDfriendTrack *friendTrack = fESDfriend->GetTrack(itr);
      if(!friendTrack) continue;
      //printf(" --- %d < %d || %p | %p -- %p \n", itr, fESDfriend->GetNumberOfTracks(), track, fESDfriend, friendTrack);
      ProcessTOF(track,friendTrack,vPos);
      if(nclusters>0)ProcessEMCAL(track,friendTrack,clusters,vPos);
      ProcessHMPID(track,friendTrack,vPos);
    }//End of track loop

  delete clusters;
  PostData(0,fOutput);
}

//________________________________________________________________________
void AliTrackComparisonESD::ProcessTOF(AliESDtrack *track, AliESDfriendTrack *friendTrack, Double_t *vPos){
  //
  // Process TPC-TOF extrapolation
  //

  //printf("Enter function!\n");
  if (track->GetTOFsignal()<=0)  return;
  if (!friendTrack) return;
  if (!friendTrack->GetTPCOut()) return;

  AliExternalTrackParam *pTPC = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
  if(!pTPC) return;

  const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
  if (!points) return;
  Int_t npoints = points->GetNPoints();
  if(npoints>1000) return; //the default value is more than 30000, why not -1???
  AliTrackPoint point;
  //
  Int_t counter=0;
  for (Int_t ipoint=0;ipoint<npoints;ipoint++){
    //cout<<"(npoints,ipoint) = ("<<npoints<<","<<ipoint<<")"<<endl;
    if(!points->GetPoint(point,ipoint)) continue;
    //Float_t xyz[3];
    //point.GetXYZ(xyz);
    //Float_t r=10;
    //Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);

    //printf("fVolumeID %d | LayerID %d\n",point.GetVolumeID(),AliGeomManager::VolUIDToLayer(point.GetVolumeID()));

    //if (r<350) continue;
    //if (r>400) continue;
    //cout<<"r="<<r<<endl;
    if(AliGeomManager::VolUIDToLayer(point.GetVolumeID())==AliGeomManager::kTOF)
      {
	counter++;
	fTOF->AddTracks(pTPC,&point,track->GetMass(),track->P(),vPos);
      }
  }
  //Printf("# of track points in TOF: %d!\n",counter);
}

//________________________________________________________________________
void AliTrackComparisonESD::ProcessEMCAL(AliESDtrack *track, AliESDfriendTrack *friendTrack, TRefArray *clusters, Double_t *vPos){
  if(clusters->GetEntriesFast()==0) return;

  Double_t rEMCal = 438;

  AliExternalTrackParam *pTPC = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
  if(!pTPC) return;

  Double_t trPos[3];
  Float_t clPos[3];
  AliExternalTrackParam *pTest = new AliExternalTrackParam(*pTPC);
  if(!AliTracker::PropagateTrackToBxByBz(pTest, rEMCal , track->GetMass(), 1 , kFALSE,0.99,-1)) return;
  if(!pTest->GetXYZ(trPos)) return;

  AliExternalTrackParam *p0=0;
  AliTrackPoint *p1=new AliTrackPoint();
  Int_t nclusters = clusters->GetEntries();
  for(Int_t icl=0; icl<nclusters; icl++)
    {
      AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(icl);
      if(!cluster) continue;
      cluster->GetPosition(clPos);
      if( TMath::Abs(clPos[0]-trPos[0])>fCutX || TMath::Abs(clPos[1]-trPos[1])>fCutY || TMath::Abs(clPos[2]-trPos[2]>fCutZ) ) continue;
      
      p0 = pTPC;
      //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]));
      p1->SetXYZ(clPos[0],clPos[1],clPos[2],0);
      //printf("Found EMCAL point!\n");
      fEMCAL->AddTracks(p0,p1,track->GetMass(),cluster->E(),vPos);
    }

  delete pTest;
  delete p1;
}

//________________________________________________________________________
void AliTrackComparisonESD::ProcessHMPID(AliESDtrack *track, AliESDfriendTrack *friendTrack, Double_t *vPos){
  //
  // Process TPC-TOF extrapolation
  //
  if (track->GetHMPIDsignal()<=0)  return;

  AliExternalTrackParam *pTPC = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
  if(!pTPC) return;

  Int_t q, nph, ch;
  Float_t x, y;
  track->GetHMPIDmip(x,y,q,nph);
  Double_t pHmp[3]={0}, pHmp3=0;
  if (track->GetOuterHmpPxPyPz(pHmp)) 
    pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);

  ch = track->GetHMPIDcluIdx()/1000000;

  AliHMPIDParam *pParam = AliHMPIDParam::Instance(); 
  TVector3 vG = pParam->Lors2Mars(ch,x,y);

  AliTrackPoint *p1 = new AliTrackPoint();
  p1->SetXYZ(vG.X(),vG.Y(),vG.Z());
  //printf("Found HMPID point!\n");
  fHMPID->AddTracks(pTPC,p1,track->GetMass(),pHmp3,vPos);
  delete p1;
}


//________________________________________________________________________
void AliTrackComparisonESD::RecalClusterPos(TRefArray *clusters, AliESDCaloCells *cells){
  Double_t iLocal[3], iGlobal[3];
  Float_t cPos[3];
  //Float_t pos[3];
  Int_t nclusters = clusters->GetEntries();
  for(Int_t icl=0; icl<nclusters; icl++)
    {
      AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(icl);
      UShort_t *absId = cluster->GetCellsAbsId();
      Int_t nCells = cluster->GetNCells();
      for(Int_t i=0;i<3;i++)cPos[i]=0;
      Double_t wTot=0;
      for(Int_t iCell=0; iCell<nCells; iCell++)
	{
	  Double_t cellEnergy = cells->GetCellAmplitude(absId[iCell]);
	  Double_t dist = 1.31*(log(cluster->E())+4.82+0.5);
	  fGeom->RelPosCellInSModule(absId[iCell],dist,iLocal[0],iLocal[1],iLocal[2]);
	  //fGeom->GetGlobal(iLocal,iGlobal,fGeom->GetSuperModuleNumber(absId[iCell]));
	  Int_t sm = fGeom->GetSuperModuleNumber(absId[iCell]);
	  //matrix[sm]->Print();
	  //cout<<"sm = "<<sm<<endl;
	  fTransMatrix[sm]->LocalToMaster(iLocal,iGlobal);

	  Double_t w = TMath::Max( 0., 4.5 + TMath::Log( cellEnergy / cluster->E() ));
	  if(w>0.0)
	    {
	      wTot += w;
	      for(Int_t i=0; i<3; i++ )
		  cPos[i] += (w*iGlobal[i]);
	    }
	}//End of cell loop   
      if(wTot>0)
	{
	  for(int i=0; i<3; i++ ) 
	    {
	      cPos[i] /= wTot;
	      cluster->SetPosition(cPos);
	    }
	}
      //cluster->GetPosition(pos);
      //printf("cluster %d pass10 pos: (%5.3f,%5.3f,%5.3f)\n",icl,pos[0],pos[1],pos[2]);
    }//End of cluster loop
}

//________________________________________________________________________
void AliTrackComparisonESD::Terminate(Option_t */*option*/) {
  //
  // Terminate
  //
  AliInfo("AliTrackComparisonESD::Terminate()\n");
  
}

//________________________________________________________________________
void AliTrackComparisonESD::FinishTaskOutput(){
  //
  // According description in AliAnalisysTask this method is call 
  // on the slaves before sending data
  //
  Terminate("slave");
  if(!fDebugOutputPath.IsNull()) { 
    RegisterDebugOutput();
  }
}

//________________________________________________________________________
Long64_t AliTrackComparisonESD::Merge(TCollection *li) {
  if(li) return 1;
  else return 1;
}

//________________________________________________________________________
void AliTrackComparisonESD::Analyze() {
  //
  // Analyze the content of the task
  //

}

//________________________________________________________________________
void AliTrackComparisonESD::RegisterDebugOutput(){
  //
  //
  //
}
Commit	Line	Data
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
	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	}