[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 "AliCalorimeterUtils.h"
#include "AliESDCaloCells.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),
   fCutR(20),
   fCaloUtil(0)
{
  //
  // 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),
   fCutR(20),
   fCaloUtil(0)
{
  //
  // 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(fEMCAL)  delete fEMCAL;  fEMCAL=0;
  if(fHMPID)  delete fHMPID;  fHMPID=0;
  if(fTOF)    delete fTOF;    fTOF=0;
  if(fOutput) fOutput->Delete();
  //if(fOutput) delete fOutput;
  if(fCaloUtil) delete fCaloUtil; fCaloUtil=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();
    }


  fCaloUtil = new AliCalorimeterUtils();
  InitCaloUtil();


  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,cells,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++){
    if(!points->GetPoint(point,ipoint)) continue;
    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);
  return;
}

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

  Double_t rEMCal = 438;
  Double_t tmp=fCutR;
  Int_t clsIndex=-1;
  TVector3 clsV,trkV;

  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;

  Double_t trPhi = TMath::ATan2(trPos[1],trPos[0])*TMath::RadToDeg();
  //printf("trPhi = %5.3f | eta = %5.3f\n",trPhi,pTest->Eta());
  if(trPhi<80 || trPhi>120) return;
  if(TMath::Abs(pTest->Eta())>0.7) 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;
      if(fCaloUtil->ClusterContainsBadChannel("EMCAL",cluster->GetCellsAbsId(),cluster->GetNCells()) ) continue;
      if(!fCaloUtil->CheckCellFiducialRegion(cluster,cells,NULL,0) ) continue;

      cluster->GetPosition(clPos);
      p0 = pTPC;
      clsV.SetXYZ(clPos[0],clPos[1],clPos[2]);
      Double_t alpha = ((int)(clsV.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
      clsV.RotateZ(-alpha);
      p0->Rotate(alpha);
      if(!AliTrackerBase::PropagateTrackToBxByBz(p0,clsV.X(), track->GetMass(), 1,kFALSE,0.99,-1)) continue;
      trkV.SetXYZ(p0->GetX(),p0->GetY(),p0->GetZ());
      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) );
      if(dist<tmp)
	{                
	  tmp=dist;
	  clsIndex=icl;
	}
    }
      
  if(clsIndex==-1) return;
  AliESDCaloCluster *cluster = (AliESDCaloCluster*) clusters->At(clsIndex);
  cluster->GetPosition(clPos);
  //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(pTPC,p1,track->GetMass(),track->P(),vPos);


  delete pTest;
  delete p1;
  return;
}

//________________________________________________________________________
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;
  return;
}


//________________________________________________________________________
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(){
  //
  //
  //
}

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