Possibility to recalculate cluster position and correct for the non linearity of...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 10 Oct 2010 18:42:28 +0000 (18:42 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 10 Oct 2010 18:42:28 +0000 (18:42 +0000)
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.h
PWG4/CaloCalib/macros/anaEMCALCalib.C
PWG4/PartCorrBase/AliCaloTrackReader.cxx
PWG4/PartCorrBase/AliCalorimeterUtils.cxx
PWG4/PartCorrBase/AliCalorimeterUtils.h

index 9f831c8..fbe10ab 100644 (file)
@@ -38,6 +38,8 @@
 #include "AliEMCALGeometry.h"
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
+#include "AliEMCALRecoUtils.h"
+
 //#include "AliEMCALAodCluster.h"
 //#include "AliEMCALCalibData.h"
 
@@ -51,9 +53,15 @@ AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(con
   fLogWeight(4.5), fCopyAOD(kFALSE), fSameSM(kFALSE), fOldAOD(kFALSE),
   fEMCALGeoName("EMCAL_FIRSTYEAR"), fNCellsFromEMCALBorder(0),
   fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(0x0),
-  fRecalibration(kFALSE),fEMCALRecalibrationFactors(), 
-  fNbins(300), fMinBin(0.), fMaxBin(300.),
-  fOutputContainer(0x0),fHmgg(0x0), fHmggDifferentSM(0x0), fhNEvents(0x0),fCuts(0x0)
+  fRecalibration(kFALSE),fEMCALRecalibrationFactors(),
+  fRecoUtils(new AliEMCALRecoUtils),
+  fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
+  fHmgg(0x0),           fHmggDifferentSM(0x0), 
+  fHOpeningAngle(0x0),  fHOpeningAngleDifferentSM(0x0),  
+  fHIncidentAngle(0x0), fHIncidentAngleDifferentSM(0x0),
+  fHAsymmetry(0x0),  fHAsymmetryDifferentSM(0x0),  
+
+  fhNEvents(0x0),fCuts(0x0)
 {
   //Named constructor which should be used.
   
@@ -66,8 +74,17 @@ AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(con
   }
   
   for(Int_t iSM=0; iSM<4; iSM++) {
-    fHmggSM[iSM]    =0;
-    fHmggPairSM[iSM]=0;
+    fHmggSM[iSM]              =0;
+    fHmggPairSM[iSM]          =0;
+    fHOpeningAngleSM[iSM]     =0;
+    fHOpeningAnglePairSM[iSM] =0;
+    fHAsymmetrySM[iSM]     =0;
+    fHAsymmetryPairSM[iSM] =0;
+    fHIncidentAngleSM[iSM]    =0;
+    fHIncidentAnglePairSM[iSM]=0;
+    fhTowerDecayPhotonHit[iSM] =0;
+    fhTowerDecayPhotonEnergy[iSM]=0;
+    fhTowerDecayPhotonAsymmetry[iSM]=0;
   }
   
   DefineOutput(1, TList::Class());
@@ -98,6 +115,9 @@ AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
                fEMCALRecalibrationFactors->Clear();
                delete  fEMCALRecalibrationFactors;
        }       
+    
+  if(fRecoUtils) delete fRecoUtils ;
+
 }
 
 //_____________________________________________________
@@ -435,7 +455,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
   //Create output container, init geometry and calibration
        
   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;  
-       
+  
   fOutputContainer = new TList();
   const Int_t buffersize = 255;
   char hname[buffersize], htitl[buffersize];
@@ -460,6 +480,37 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
   fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
   fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
   fOutputContainer->Add(fHmggDifferentSM);
+
+  fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
+  fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
+  fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+  fOutputContainer->Add(fHOpeningAngle);
+  
+  fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
+  fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
+  fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+  fOutputContainer->Add(fHOpeningAngleDifferentSM);
+  
+  fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
+  fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
+  fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fHIncidentAngle);
+  
+  fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
+  fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
+  fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fHIncidentAngleDifferentSM);
+  
+  fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
+  fHAsymmetry->SetXTitle("a");
+  fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+  fOutputContainer->Add(fHAsymmetry);
+  
+  fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
+  fHAsymmetryDifferentSM->SetXTitle("a");
+  fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+  fOutputContainer->Add(fHAsymmetryDifferentSM);
+  
   
   TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
   
@@ -470,19 +521,84 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
     fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
     fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
     fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
-    
     fOutputContainer->Add(fHmggSM[iSM]);
     
-    
     snprintf(hname,buffersize, "hmgg_PairSM%d",iSM);
     snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair: %s",pairname[iSM].Data());
     fHmggPairSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
     fHmggPairSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
     fHmggPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
-    
     fOutputContainer->Add(fHmggPairSM[iSM]);
+    
+    
+    snprintf(hname, buffersize, "hopang_SM%d",iSM);
+    snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
+    fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
+    fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
+    fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+    fOutputContainer->Add(fHOpeningAngleSM[iSM]);
+    
+    snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
+    snprintf(htitl,buffersize, "Opening angle for SM pair: %s",pairname[iSM].Data());
+    fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
+    fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
+    fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+    fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);    
+    
+    snprintf(hname, buffersize, "hinang_SM%d",iSM);
+    snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
+    fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
+    fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
+    fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fHIncidentAngleSM[iSM]);
+    
+    snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
+    snprintf(htitl,buffersize, "Incident angle for SM pair: %s",pairname[iSM].Data());
+    fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
+    fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
+    fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);   
+    
+    snprintf(hname, buffersize, "hasym_SM%d",iSM);
+    snprintf(htitl, buffersize, "asymmetry for super mod %d",iSM);
+    fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
+    fHAsymmetrySM[iSM]->SetXTitle("a");
+    fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+    fOutputContainer->Add(fHAsymmetrySM[iSM]);
+    
+    snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
+    snprintf(htitl,buffersize, "Asymmetry for SM pair: %s",pairname[iSM].Data());
+    fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
+    fHAsymmetryPairSM[iSM]->SetXTitle("a");
+    fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+    fOutputContainer->Add(fHAsymmetryPairSM[iSM]);    
+    
+    
+    Int_t colmax = 48;
+    Int_t rowmax = 24;
+    
+    fhTowerDecayPhotonHit[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),Form("Entries in grid of cells in Module %d",iSM), 
+                                      colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
+    fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
+    fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
+    fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
+    
+    fhTowerDecayPhotonEnergy[iSM]  = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),Form("Accumulated energy in grid of cells in Module %d",iSM), 
+                                       colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
+    fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
+    fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
+    fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
+    
+    fhTowerDecayPhotonAsymmetry[iSM]  = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),Form("Accumulated asymmetry in grid of cells in Module %d",iSM), 
+                                               colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
+    fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
+    fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
+    fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
+    
   }  
   
+  
+  
   fhNEvents        = new TH1I("hNEvents", "Number of analyzed events"   , 1 , 0 , 1  ) ;
   fOutputContainer->Add(fhNEvents);
 
@@ -574,6 +690,9 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
        
   TLorentzVector p1;
   TLorentzVector p2;
+//  TLorentzVector p11;
+//  TLorentzVector p22;
+
   TLorentzVector p12;
   
   TRefArray * caloClustersArr  = new TRefArray();
@@ -585,7 +704,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
   
   // EMCAL cells
   AliAODCaloCells *emCells = aod->GetEMCALCells();
-  
+   
   // loop over EMCAL clusters
   for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
     
@@ -629,12 +748,25 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
     }
     
     //clu1.GetMomentum(p1,v);
-    c1->GetMomentum(p1,v);
-    
-    //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
-    //MaxEnergyCellPos(emCells,&clu1,iSupMod1,ieta1,iphi1);
-    MaxEnergyCellPos(emCells,c1,iSupMod1,ieta1,iphi1);
     
+    // Correct Non-Linearity
+    c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
+    //printf("\t  e1 org %2.3f, e1 cor  %2.3f \n",e1ii,c1->E());
+
+    //c1->GetMomentum(p1,v);
+    //printf("\t cor: e %2.3f, pt %2.3f, px %2.3f, py %2.3f, pz %2.3f\n",p1.E(), p1.Pt(),p1.Px(),p1.Py(),p1.Pz());
+
+    //Get tower with maximum energy and fill in the end the pi0 histogram for this cell, recalculate cluster position and recalibrate    
+    //Float_t pos[3];
+    //c1->GetPosition(pos);
+    //printf("Before Alignment: e %2.4f, x %2.4f, y %2.4f , z %2.4f\n",c1->E(),pos[0], pos[1],pos[2]);
+    GetMaxEnergyCellPosAndClusterPos(emCells,c1,iSupMod1,ieta1,iphi1);
+    //printf("i1 %d, corr1 %2.3f, e1 %2.3f, , ecorr1 %2.3f, pt %2.3f, px %2.3f, py %2.3f, pz %2.3f,\n",iClu, 1./corrFac, e1, p1.E(), p1.Pt(),p1.Px(),p1.Py(),p1.Pz());
+    //c1->GetPosition(pos);
+    //printf("After Alignment: e %2.4f, x %2.4f, y %2.4f , z %2.4f\n",c1->E(),pos[0], pos[1],pos[2]);
+
+    c1->GetMomentum(p1,v);
+
     for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
       if(c2->IsEqual(c1)) continue;
@@ -657,16 +789,20 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
       
       Float_t e2ii = c2->E();
       
-      //clu2.GetMomentum(p2,v);
+      //Correct Non-Linearity
+      c2->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c2));
+
+      //Get tower with maximum energy and fill in the end the pi0 histogram for this cell, recalculate cluster position and recalibrate    
+      GetMaxEnergyCellPosAndClusterPos(emCells,c2,iSupMod2,ieta2,iphi2);
+  
       c2->GetMomentum(p2,v);
-      
-      //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
-      //MaxEnergyCellPos(emCells,&clu2,iSupMod2,ieta2,iphi2);
-      MaxEnergyCellPos(emCells,c2,iSupMod2,ieta2,iphi2);
-      
+
       p12 = p1+p2;
       Float_t invmass = p12.M()*1000; 
-      
+      //printf("*** mass %f\n",invmass);
+      Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
+      //printf("asymmetry %f\n",asym);
+
       if(invmass < fMaxBin && invmass > fMinBin){
         
         //Check if cluster is in fidutial region, not too close to borders
@@ -684,8 +820,80 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
           if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) fHmggPairSM[1]->Fill(invmass,p12.Pt()); 
           if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) fHmggPairSM[2]->Fill(invmass,p12.Pt()); 
           if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) fHmggPairSM[3]->Fill(invmass,p12.Pt()); 
+          
+          if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
+            
+            //Opening angle of 2 photons
+            Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
+            //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
+
+            //Inciden angle of each photon
+            //Float_t * posSM1cen = RecalculatePosition(11.5, 23.5, p1.E(),0, iSupMod1); 
+            //Float_t * posSM2cen = RecalculatePosition(11.5, 23.5, p2.E(),0, iSupMod2); 
+            Float_t posSM1cen[3]={0.,0.,0.};
+            fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, p1.E(),iSupMod1,0,fRecoUtils->GetMisalShiftArray(),posSM1cen); 
+            Float_t posSM2cen[3]={0.,0.,0.}; 
+            fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, p2.E(),iSupMod2,0,fRecoUtils->GetMisalShiftArray(),posSM2cen); 
+            //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
+            //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
+            
+            TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]); 
+            TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]); 
+            Float_t inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
+            Float_t inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
+            //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
+            
+            fHOpeningAngle ->Fill(opangle,p12.Pt()); 
+            fHIncidentAngle->Fill(inangle1,p1.Pt()); 
+            fHIncidentAngle->Fill(inangle2,p2.Pt()); 
+            fHAsymmetry    ->Fill(asym,p12.Pt()); 
+
+            if(iSupMod1==iSupMod2) {
+              fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
+              fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
+              fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
+              fHAsymmetrySM[iSupMod1]    ->Fill(asym,p12.Pt());
+            }
+            else{      
+              fHOpeningAngleDifferentSM  ->Fill(opangle,p12.Pt());
+              fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
+              fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
+              fHAsymmetryDifferentSM     ->Fill(asym,p12.Pt());
+            }
+            
+            if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
+              fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt()); 
+              fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
+              fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
+              fHAsymmetryPairSM[0]    ->Fill(asym,p12.Pt());
+
+            } 
+            if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
+              fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt()); 
+              fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
+              fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
+              fHAsymmetryPairSM[1]    ->Fill(asym,p12.Pt());
+
+            }
+            
+            if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
+              fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt()); 
+              fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
+              fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
+              fHAsymmetryPairSM[2]    ->Fill(asym,p12.Pt());
+
+
+            }
+            if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
+              fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt()); 
+              fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
+              fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
+              fHAsymmetryPairSM[3]    ->Fill(asym,p12.Pt());
+            }
+              
+          }// pair in 100 < mass < 160
         
-        }
+        }//in acceptance cuts
         
         //In case of filling only channels with second cluster in same SM
         if(fSameSM && iSupMod1!=iSupMod2) continue;
@@ -693,6 +901,17 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
         if (fGroupNCells == 0){
             fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
             fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
+          
+            if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
+              fhTowerDecayPhotonHit      [iSupMod1]->Fill(ieta1,iphi1);
+              fhTowerDecayPhotonEnergy   [iSupMod1]->Fill(ieta1,iphi1,p1.E());
+              fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
+              
+              fhTowerDecayPhotonHit      [iSupMod2]->Fill(ieta2,iphi2);
+              fhTowerDecayPhotonEnergy   [iSupMod2]->Fill(ieta2,iphi2,p2.E());
+              fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
+
+            }// pair in mass of pi0
         }      
         else  {
           //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
@@ -795,7 +1014,7 @@ Bool_t AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion(AliVCluste
 }      
 
 //__________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
+void AliAnalysisTaskEMCALPi0CalibSelection::GetMaxEnergyCellPosAndClusterPos(AliVCaloCells* cells, AliVCluster* clu, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
 {
   //For a given CaloCluster calculates the absId of the cell 
   //with maximum energy deposit.
@@ -811,18 +1030,33 @@ void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* co
   Int_t iIphi   = -1;
   Int_t iIeta   = -1;
        
+  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
+  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
+  Bool_t  areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
+  Int_t   startingSM = -1;
+  
   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
     cellAbsId = clu->GetCellAbsId(iDig);
     fraction  = clu->GetCellAmplitudeFraction(iDig);
     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+    Int_t imodrc   = -1, iphirc  = -1, ietarc  =-1;
+    Int_t iTowerrc = -1, iIphirc = -1, iIetarc =-1;
+    fEMCALGeo->GetCellIndex(cellAbsId,imodrc,iTowerrc,iIphirc,iIetarc); 
+    fEMCALGeo->GetCellPhiEtaIndexInSModule(imodrc,iTowerrc,iIphirc, iIetarc,iphirc,ietarc);                    
+    if     (iDig==0)  startingSM = imodrc;
+    else if(imodrc != startingSM) areInSameSM = kFALSE;
+
     if(IsRecalibrationOn()) {
-      Int_t imodrc   = -1, iphirc  = -1, ietarc  =-1;
-      Int_t iTowerrc = -1, iIphirc = -1, iIetarc =-1;
-      fEMCALGeo->GetCellIndex(cellAbsId,imodrc,iTowerrc,iIphirc,iIetarc); 
-      fEMCALGeo->GetCellPhiEtaIndexInSModule(imodrc,iTowerrc,iIphirc, iIetarc,iphirc,ietarc);                  
       recalFactor = GetEMCALChannelRecalibrationFactor(imodrc,ietarc,iphirc);
     }
     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
+    
+    weight = TMath::Log(eCell/clEnergy) + 4;
+    if(weight < 0) weight = 0;
+    totalWeight += weight;
+    weightedCol += ietarc*weight;
+    weightedRow += iphirc*weight;
+    
     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
     
     if(eCell > eMax)  { 
@@ -830,14 +1064,34 @@ void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* co
       maxId = cellAbsId;
       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
     }
-  }
+  }// cell loop
   
   //Get from the absid the supermodule, tower and eta/phi numbers
   fEMCALGeo->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta); 
   //Gives SuperModule and Tower numbers
   fEMCALGeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                        iIphi, iIeta,iphi,ieta);    
+                                        iIphi, iIeta,iphi,ieta); 
   
+  Float_t xyzNew[3];
+  if(areInSameSM == kTRUE) {
+    //printf("In Same SM\n");
+    weightedCol = weightedCol/totalWeight;
+    weightedRow = weightedRow/totalWeight;
+    
+    //Float_t *xyzNew = RecalculatePosition(weightedRow, weightedCol, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
+    fEMCALGeo->RecalculateTowerPosition(weightedRow, weightedCol, iSupMod, clEnergy, 0, //1 = electrons, 0 photons
+                                        fRecoUtils->GetMisalShiftArray(), xyzNew);
+  }
+  else {
+    //printf("In Different SM\n");
+    //Float_t *xyzNew = RecalculatePosition(iphi,        ieta,        clEnergy, 0, iSupMod); //1 = electrons, 0 photons
+    fEMCALGeo->RecalculateTowerPosition(iphi, ieta, iSupMod, clEnergy, 0, //1 = electrons, 0 photons
+                                        fRecoUtils->GetMisalShiftArray(), xyzNew);
+    
+  }
+
+  clu->SetPosition(xyzNew);
+
   //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
   
 }
index f7ddf98..ecbf20b 100644 (file)
@@ -21,16 +21,22 @@ class AliAODCaloCluster;
 class AliAODCaloCells;
 //class AliEMCALCalibData ;
 #include "AliEMCALGeoParams.h"
+class AliEMCALRecoUtils;
 
 class AliAnalysisTaskEMCALPi0CalibSelection : public AliAnalysisTaskSE
 {
 public:
 
   AliAnalysisTaskEMCALPi0CalibSelection(const char* name);
-  AliAnalysisTaskEMCALPi0CalibSelection(const AliAnalysisTaskEMCALPi0CalibSelection&); 
-  AliAnalysisTaskEMCALPi0CalibSelection& operator=(const AliAnalysisTaskEMCALPi0CalibSelection&); 
   virtual ~AliAnalysisTaskEMCALPi0CalibSelection();
 
+private:
+  
+  AliAnalysisTaskEMCALPi0CalibSelection(const AliAnalysisTaskEMCALPi0CalibSelection&); 
+  AliAnalysisTaskEMCALPi0CalibSelection& operator=(const AliAnalysisTaskEMCALPi0CalibSelection&); 
+  
+public:
+  
   // Implementation of interface methods
   virtual void UserCreateOutputObjects();
   virtual void UserExec(Option_t * opt);
@@ -92,7 +98,7 @@ public:
   void SwitchOffRecalibration()   {fRecalibration = kFALSE ; }
        
   void InitEMCALRecalibrationFactors() ;
-       
+  
   Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const { 
        if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow); 
        else return 1;}
@@ -108,12 +114,14 @@ public:
   void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecalibrationFactors = map;}
   Float_t RecalibrateClusterEnergy(AliAODCaloCluster* cluster, AliAODCaloCells * cells);
        
+  void SetEMCALRecoUtils(AliEMCALRecoUtils * ru) {fRecoUtils = ru;}
+  AliEMCALRecoUtils* GetEMCALRecoUtils() const {return fRecoUtils;}
+  
   void SetInvariantMassHistoBinRange(Int_t nBins, Float_t minbin, Float_t maxbin){
        fNbins = nBins; fMinBin = minbin; fMaxBin = maxbin; }
-       
+         
 private:
-
-  void MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& iSM, Int_t& ieta, Int_t& iphi);
+  void GetMaxEnergyCellPosAndClusterPos(AliVCaloCells* cells, AliVCluster* clu, Int_t& iSM, Int_t& ieta, Int_t& iphi);
 
 private:
 
@@ -136,6 +144,8 @@ private:
   Bool_t     fRecalibration;             // Switch on or off the recalibration
   TObjArray *fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL                 
  
+  AliEMCALRecoUtils * fRecoUtils;  // Access to reconstruction utilities
+  
   //Output histograms  
   Int_t   fNbins;  // N       mass bins of invariant mass histograms
   Float_t fMinBin; // Minimum mass bins of invariant mass histograms
@@ -143,14 +153,35 @@ private:
 
   TList*  fOutputContainer; //!histogram container
   TH1F*   fHmpi0[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];//! two-cluster inv. mass assigned to each cell.
-  TH2F*   fHmgg;            //! two-cluster inv.mass vt pt of pair
-  TH2F*   fHmggDifferentSM; //! two-cluster inv.mass vt pt of pair, each cluster in different SM
+
+  TH2F*   fHmgg;            //! two-cluster inv.mass vs pt of pair
+  TH2F*   fHmggDifferentSM; //! two-cluster inv.mass vs pt of pair, each cluster in different SM
   TH2F*   fHmggSM[4];       //! two-cluster inv.mass per SM
   TH2F*   fHmggPairSM[4];   //! two-cluster inv.mass per Pair
+  
+  TH2F*   fHOpeningAngle;            //! two-cluster opening angle vs pt of pair, with mass close to pi0
+  TH2F*   fHOpeningAngleDifferentSM; //! two-cluster opening angle vs pt of pair, each cluster in different SM, with mass close to pi0
+  TH2F*   fHOpeningAngleSM[4];       //! two-cluster opening angle vs pt per SM,with mass close to pi0
+  TH2F*   fHOpeningAnglePairSM[4];   //! two-cluster opening angle vs pt per Pair,with mass close to pi0
+
+  TH2F*   fHIncidentAngle;            //! cluster incident angle vs pt of pair, with mass close to pi0
+  TH2F*   fHIncidentAngleDifferentSM; //! cluster incident angle vs pt of pair, each cluster in different SM, with mass close to pi0
+  TH2F*   fHIncidentAngleSM[4];       //! cluster incident angle vs pt per SM,with mass close to pi0
+  TH2F*   fHIncidentAnglePairSM[4];   //! cluster incident angle vs pt per Pair,with mass close to pi0
+  
+  TH2F*   fHAsymmetry;            //! two-cluster asymmetry vs pt of pair, with mass close to pi0
+  TH2F*   fHAsymmetryDifferentSM; //! two-cluster asymmetry vs pt of pair, each cluster in different SM, with mass close to pi0
+  TH2F*   fHAsymmetrySM[4];       //! two-cluster asymmetry vs pt per SM,with mass close to pi0
+  TH2F*   fHAsymmetryPairSM[4];   //! two-cluster asymmetry vs pt per Pair,with mass close to pi0
+  
+  TH2F*   fhTowerDecayPhotonHit[4] ;       //! Cells ordered in column/row for different module, number of times a decay photon hits
+  TH2F*   fhTowerDecayPhotonEnergy[4] ;    //! Cells ordered in column/row for different module, accumulated energy in the tower by decay photons.
+  TH2F*   fhTowerDecayPhotonAsymmetry[4] ; //! Cells ordered in column/row for different module, accumulated asymmetry in the tower by decay photons.
+
   TH1I*   fhNEvents;        //! Number of events counter histogram
   TList * fCuts ;           //! List with analysis cuts
 
-  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,5);
+  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,7);
 
 };
 
index 82dbf3d..e1e91b9 100644 (file)
@@ -28,13 +28,15 @@ char * kXML = "collection.xml";
 //---------------------------------------------------------------------------
 
 const TString kInputData = "AOD"; //ESD, AOD, MC
-TString kTreeName = "aodTree";
+TString kTreeName = "esdTree";
 Bool_t copy = kFALSE;
 
 void anaEMCALCalib(Int_t mode=mLocal)
 {
   // Main
-
+  char cmd[200] ; 
+  sprintf(cmd, ".! rm -rf aod.root pi0calib.root") ; 
+  gROOT->ProcessLine(cmd) ; 
   //--------------------------------------------------------------------
   // Load analysis libraries
   // Look at the method below, 
@@ -100,16 +102,21 @@ void anaEMCALCalib(Int_t mode=mLocal)
     //-------------------------------------------------------------------------
     AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
-    
-    gROOT->LoadMacro("AddTaskPhysicsSelection.C");
-    AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
+
     
     // ESD filter task
-    if(kInputData == "ESD" && !copy){
-      gROOT->LoadMacro("AddTaskESDFilter.C");
-      AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
+    if(kInputData == "ESD"){
+
+      gROOT->LoadMacro("AddTaskPhysicsSelection.C");
+      AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
+      if(!copy){
+       gROOT->LoadMacro("AddTaskESDFilter.C");
+       AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
+      }
     }
     
+
+   
     AliAnalysisTaskEMCALPi0CalibSelection * pi0calib = new AliAnalysisTaskEMCALPi0CalibSelection ("EMCALPi0Calibration");
     //pi0calib->SetDebugLevel(10); 
     pi0calib->CopyAOD(copy);
@@ -119,7 +126,27 @@ void anaEMCALCalib(Int_t mode=mLocal)
     pi0calib->SetNCellsGroup(0);
     pi0calib->SwitchOnBadChannelsRemoval();
     pi0calib->SwitchOffSameSM();
-    pi0calib->SwitchOffOldAODs();
+    pi0calib->SwitchOnOldAODs();
+    pi0calib->SetNumberOfCellsFromEMCALBorder(1);
+    AliEMCALRecoUtils * reco = pi0calib->GetEMCALRecoUtils();
+    reco->SetMisalShift(0,0.8);  reco->SetMisalShift(1,8.3); reco->SetMisalShift(2,1.);
+    reco->SetMisalShift(3,-7.5); reco->SetMisalShift(4,7.5); reco->SetMisalShift(5,2.);
+    reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaGamma);
+    reco->SetNonLinearityParam(0,1.04);     reco->SetNonLinearityParam(1,-0.1445);
+    reco->SetNonLinearityParam(2,1.046);    
+
+//     reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaConversion);
+//     reco->SetNonLinearityParam(0,1.033);     reco->SetNonLinearityParam(1,0.0566186);
+//     reco->SetNonLinearityParam(2,0.982133);    
+
+
+//      reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MC);
+//      reco->SetNonLinearityParam(0,1.001);     reco->SetNonLinearityParam(1,-0.01264);
+//      reco->SetNonLinearityParam(2,-0.03632);    
+//      reco->SetNonLinearityParam(3,0.1798);     reco->SetNonLinearityParam(4,-0.522);
+
+    reco->Print("");
+
     // SM0
     pi0calib->SetEMCALChannelStatus(0,3,13);  pi0calib->SetEMCALChannelStatus(0,44,1); pi0calib->SetEMCALChannelStatus(0,3,13); 
     pi0calib->SetEMCALChannelStatus(0,20,7);  pi0calib->SetEMCALChannelStatus(0,38,2);   
@@ -140,6 +167,17 @@ void anaEMCALCalib(Int_t mode=mLocal)
 
     mgr->AddTask(pi0calib);
 
+   pi0calib->SwitchOnRecalibration();
+  TFile * f = new TFile("RecalibrationFactors.root","read");
+  TH2F * h0 = (TH2F*)f->Get("EMCALRecalFactors_SM0")->Clone();
+  TH2F * h1 = (TH2F*)f->Get("EMCALRecalFactors_SM1")->Clone();
+  TH2F * h2 = (TH2F*)f->Get("EMCALRecalFactors_SM2")->Clone();
+  TH2F * h3 = (TH2F*)f->Get("EMCALRecalFactors_SM3")->Clone();
+
+  pi0calib->SetEMCALChannelRecalibrationFactors(0,h0);
+  pi0calib->SetEMCALChannelRecalibrationFactors(1,h1);
+  pi0calib->SetEMCALChannelRecalibrationFactors(2,h2);
+  pi0calib->SetEMCALChannelRecalibrationFactors(3,h3);
     
     // Create containers for input/output
     AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
@@ -148,7 +186,7 @@ void anaEMCALCalib(Int_t mode=mLocal)
       mgr->CreateContainer("pi0calib", TList::Class(), AliAnalysisManager::kOutputContainer, "pi0calib.root");
     
     AliAnalysisDataContainer *cout_cuts = mgr->CreateContainer("Cuts", TList::Class(), 
-                                                              AliAnalysisManager::kParamContainer, "pi0calib.root");
+                                                              AliAnalysisManager::kOutputContainer, "pi0calib.root");
     
     mgr->ConnectInput  (pi0calib,     0, cinput1);
     if(copy) mgr->ConnectOutput (pi0calib, 0, coutput1 );
@@ -202,6 +240,7 @@ void  LoadLibraries(const anaModes mode) {
     gSystem->Load("libANALYSIS.so");
     gSystem->Load("libANALYSISalice.so");
     gSystem->Load("libANALYSISalice.so");
+    TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
     gSystem->Load("libPWG4CaloCalib.so");
     
     /*
index d0a76a5..85a5220 100755 (executable)
@@ -43,6 +43,7 @@
 #include "AliVParticle.h"
 #include "AliMixedEvent.h"
 #include "AliESDtrack.h"
+#include "AliEMCALRecoUtils.h"
 
 ClassImp(AliCaloTrackReader)
   
@@ -66,104 +67,12 @@ ClassImp(AliCaloTrackReader)
     fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
     fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0), 
     fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL), 
-    fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE)
-{
+    fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE){
   //Ctor
   
   //Initialize parameters
   InitParameters();
 }
-/*
-//____________________________________________________________________________
-AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & reader) :   
-  TObject(reader), fEventNumber(reader.fEventNumber), fCurrentFileName(reader.fCurrentFileName), 
-  fDataType(reader.fDataType), fDebug(reader.fDebug),
-  fFiducialCut(reader.fFiducialCut),
-  fComparePtHardAndJetPt(reader.fComparePtHardAndJetPt),
-  fPtHardAndJetPtFactor(reader.fPtHardAndJetPtFactor),
-  fCTSPtMin(reader.fCTSPtMin), fEMCALPtMin(reader.fEMCALPtMin),fPHOSPtMin(reader.fPHOSPtMin), 
-  fAODBranchList(new TList()),
-  fAODCTS(new TObjArray(*reader.fAODCTS)),  
-  fAODEMCAL(new TObjArray(*reader.fAODEMCAL)),
-  fAODPHOS(new TObjArray(*reader.fAODPHOS)),
-  fEMCALCells(new TNamed(*reader.fEMCALCells)),
-  fPHOSCells(new TNamed(*reader.fPHOSCells)),
-  fInputEvent(reader.fInputEvent), fOutputEvent(reader.fOutputEvent), fMC(reader.fMC),
-  fFillCTS(reader.fFillCTS),fFillEMCAL(reader.fFillEMCAL),fFillPHOS(reader.fFillPHOS),
-  fFillEMCALCells(reader.fFillEMCALCells),fFillPHOSCells(reader.fFillPHOSCells),
-  fSecondInputAODTree(reader.fSecondInputAODTree), 
-  fSecondInputAODEvent(reader.fSecondInputAODEvent),
-  fSecondInputFileName(reader.fSecondInputFileName), 
-  fSecondInputFirstEvent(reader.fSecondInputFirstEvent),
-  fAODCTSNormalInputEntries(reader.fAODCTSNormalInputEntries), 
-  fAODEMCALNormalInputEntries(reader.fAODEMCALNormalInputEntries), 
-  fAODPHOSNormalInputEntries(reader.fAODPHOSNormalInputEntries),
-  fTrackStatus(reader.fTrackStatus),
-  fReadStack(reader.fReadStack), fReadAODMCParticles(reader.fReadAODMCParticles),
-  fFiredTriggerClassName(reader.fFiredTriggerClassName),
-  fAnaLED(reader.fAnaLED), 
-  fTaskName(reader.fTaskName),
-  fCaloUtils(new AliCalorimeterUtils(*reader.fCaloUtils))
-{
-  // cpy ctor  
-}
-*/
-//_________________________________________________________________________
-//AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & source)
-//{
-//  // assignment operator
-//  
-//  if(&source == this) return *this;
-//  
-//  fDataType    = source.fDataType ;
-//  fDebug       = source.fDebug ;
-//  fEventNumber = source.fEventNumber ;
-//  fCurrentFileName = source.fCurrentFileName ;
-//  fFiducialCut = source.fFiducialCut;
-//     
-//  fComparePtHardAndJetPt = source.fComparePtHardAndJetPt;
-//  fPtHardAndJetPtFactor  = source.fPtHardAndJetPtFactor;
-//     
-//  fCTSPtMin    = source.fCTSPtMin ;
-//  fEMCALPtMin  = source.fEMCALPtMin ;
-//  fPHOSPtMin   = source.fPHOSPtMin ; 
-//  
-//  fAODCTS     = new TObjArray(*source.fAODCTS) ;
-//  fAODEMCAL   = new TObjArray(*source.fAODEMCAL) ;
-//  fAODPHOS    = new TObjArray(*source.fAODPHOS) ;
-//  fEMCALCells = new TNamed(*source.fEMCALCells) ;
-//  fPHOSCells  = new TNamed(*source.fPHOSCells) ;
-//
-//  fInputEvent  = source.fInputEvent;
-//  fOutputEvent = source.fOutputEvent;
-//  fMC          = source.fMC;
-//  
-//  fFillCTS        = source.fFillCTS;
-//  fFillEMCAL      = source.fFillEMCAL;
-//  fFillPHOS       = source.fFillPHOS;
-//  fFillEMCALCells = source.fFillEMCALCells;
-//  fFillPHOSCells  = source.fFillPHOSCells;
-//
-//  fSecondInputAODTree    = source.fSecondInputAODTree;
-//  fSecondInputAODEvent   = source.fSecondInputAODEvent;
-//  fSecondInputFileName   = source.fSecondInputFileName;
-//  fSecondInputFirstEvent = source.fSecondInputFirstEvent;
-//
-//  fAODCTSNormalInputEntries   = source.fAODCTSNormalInputEntries; 
-//  fAODEMCALNormalInputEntries = source.fAODEMCALNormalInputEntries; 
-//  fAODPHOSNormalInputEntries  = source.fAODPHOSNormalInputEntries;
-//     
-//  fTrackStatus        = source.fTrackStatus;
-//  fReadStack          = source.fReadStack;
-//  fReadAODMCParticles = source.fReadAODMCParticles;  
-//     
-//  fDeltaAODFileName   = source.fDeltaAODFileName;
-//     
-//  fFiredTriggerClassName = source.fFiredTriggerClassName  ;
-//     
-//  return *this;
-//  
-//}
 
 //_________________________________
 AliCaloTrackReader::~AliCaloTrackReader() {
@@ -739,10 +648,28 @@ void AliCaloTrackReader::FillInputEMCAL() {
             printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
           
+          Float_t pos[3];
+          clus->GetPosition(pos);
+          //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
+          
           //Recalibrate the cluster energy 
           if(GetCaloUtils()->IsRecalibrationOn()) {
-            Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetEMCALCells());
+            Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
             clus->SetE(energy);
+            //printf("Recalibrated Energy %f\n",clus->E());  
+          }
+          
+          //Correct non linearity
+          if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn()){
+            GetCaloUtils()->CorrectClusterEnergy(clus) ;
+            //printf("Linearity Corrected Energy %f\n",clus->E());  
+          }
+          
+          //Recalculate cluster position
+          if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
+            GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus); 
+            clus->GetPosition(pos);
+            //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
           }
           
           if (fMixedEvent) {
index 1a7c6d1..588d08d 100755 (executable)
@@ -35,6 +35,7 @@
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
 #include "AliMixedEvent.h"
+#include "AliEMCALRecoUtils.h"
 
 
 ClassImp(AliCalorimeterUtils)
@@ -50,7 +51,9 @@ ClassImp(AliCalorimeterUtils)
     fEMCALBadChannelMap(0x0),fPHOSBadChannelMap(0x0), 
     fNCellsFromEMCALBorder(0), fNCellsFromPHOSBorder(0), 
     fNoEMCALBorderAtEta0(kFALSE),fRecalibration(kFALSE),
-    fEMCALRecalibrationFactors(), fPHOSRecalibrationFactors()
+    fEMCALRecalibrationFactors(), fPHOSRecalibrationFactors(),
+    fEMCALRecoUtils(new AliEMCALRecoUtils),fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE)
+
 {
   //Ctor
   
@@ -83,6 +86,8 @@ AliCalorimeterUtils::~AliCalorimeterUtils() {
     delete  fPHOSRecalibrationFactors;
   }
        
+  if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
+  
 }
 
 //_______________________________________________________________
@@ -266,6 +271,12 @@ Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort
 }
 
 //____________________________________________________________________________________________________________________________________________________
+void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus){
+  // Correct cluster energy non linearity
+  clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
+}
+
+//____________________________________________________________________________________________________________________________________________________
 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
 {
        //Get the EMCAL/PHOS module number that corresponds to this particle
@@ -567,6 +578,8 @@ void AliCalorimeterUtils::Print(const Option_t * opt) const
         fNCellsFromEMCALBorder, fNCellsFromPHOSBorder);
   if(fNoEMCALBorderAtEta0) printf("Do not remove EMCAL clusters at Eta = 0\n");
   printf("Recalibrate Clusters? %d\n",fRecalibration);
+  printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
+  printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
 
   printf("    \n") ;
 } 
@@ -661,3 +674,81 @@ void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEven
 
 }
 
+//________________________________________________________________
+void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu){
+  
+  //Recalculate EMCAL cluster position
+  
+  Double_t eMax       = -1.;
+  Double_t eCell      = -1.;
+  Float_t  fraction   = 1.;
+  Int_t    cellAbsId  = -1;
+  Float_t recalFactor = 1.;
+       
+  Int_t maxId   = -1;
+  Int_t imod   = -1, iphi  = -1, ieta  =-1;
+  Int_t iTower = -1, iIphi = -1, iIeta =-1;
+
+  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
+  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
+  Bool_t  areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
+  Int_t   startingSM = -1;
+  
+  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
+    cellAbsId = clu->GetCellAbsId(iDig);
+    fraction  = clu->GetCellAmplitudeFraction(iDig);
+    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+       fEMCALGeo->GetCellIndex(cellAbsId,imod,iTower,iIphi,iIeta); 
+    fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);                        
+    if     (iDig==0)  startingSM = imod;
+    else if(imod != startingSM) areInSameSM = kFALSE;
+    
+    if(IsRecalibrationOn()) {
+      recalFactor = GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
+    }
+    eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
+    
+    weight = TMath::Log(eCell/clEnergy) + 4;
+    if(weight < 0) weight = 0;
+    totalWeight += weight;
+    weightedCol += ieta*weight;
+    weightedRow += iphi*weight;
+    
+    //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
+    
+    if(eCell > eMax)  { 
+      eMax  = eCell; 
+      maxId = cellAbsId;
+      //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
+    }
+  }// cell loop
+  
+  //Get from the absid the supermodule, tower and eta/phi numbers
+  fEMCALGeo->GetCellIndex(maxId,imod,iTower,iIphi,iIeta); 
+  //Gives SuperModule and Tower numbers
+  fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,
+                                         iIphi, iIeta,iphi,ieta); 
+  
+  Float_t xyzNew[3];
+  if(areInSameSM == kTRUE) {
+    //printf("In Same SM\n");
+    weightedCol = weightedCol/totalWeight;
+    weightedRow = weightedRow/totalWeight;
+    
+    //Float_t *xyzNew = RecalculatePosition(weightedRow, weightedCol, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
+    fEMCALGeo->RecalculateTowerPosition(weightedRow, weightedCol, imod, clEnergy, 0, //1 = electrons, 0 photons
+                                        fEMCALRecoUtils->GetMisalShiftArray(), xyzNew);
+  }
+  else {
+    //printf("In Different SM\n");
+    //Float_t *xyzNew = RecalculatePosition(iphi,        ieta,        clEnergy, 0, iSupMod); //1 = electrons, 0 photons
+    fEMCALGeo->RecalculateTowerPosition(iphi, ieta, imod, clEnergy, 0, //1 = electrons, 0 photons
+                                        fEMCALRecoUtils->GetMisalShiftArray(), xyzNew);
+    
+  }
+  
+  clu->SetPosition(xyzNew);
+  
+  //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
+  
+}
index 2d06ddc..2a08022 100755 (executable)
@@ -26,6 +26,7 @@ class AliVCluster;
 class AliVCaloCells;
 #include "AliPHOSGeoUtils.h"
 #include "AliEMCALGeoUtils.h"
+class AliEMCALRecoUtils;
 
 class AliCalorimeterUtils : public TObject {
 
@@ -144,6 +145,19 @@ class AliCalorimeterUtils : public TObject {
 
   Float_t RecalibrateClusterEnergy(AliVCluster* cluster, AliVCaloCells * cells);
 
+  void SetEMCALRecoUtils(AliEMCALRecoUtils * ru) {fEMCALRecoUtils = ru;}
+  AliEMCALRecoUtils* GetEMCALRecoUtils() const {return fEMCALRecoUtils;}
+  
+  Bool_t IsCorrectionOfClusterEnergyOn()  const    { return fCorrectELinearity ; }
+  void SwitchOnCorrectClusterLinearity()         { fCorrectELinearity = kTRUE; } 
+  void SwitchOffCorrectClusterLinearity()        { fCorrectELinearity = kFALSE; } 
+  void CorrectClusterEnergy(AliVCluster *cl);
+  
+  Bool_t IsRecalculationOfClusterPositionOn()  const { return fRecalculatePosition ; }
+  void SwitchOnRecalculateClusterPosition()      { fRecalculatePosition = kTRUE; } 
+  void SwitchOffRecalculateClusterPosition()     { fRecalculatePosition = kFALSE; } 
+  void RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu);
+
  private:
 
   Int_t              fDebug;                 //  Debugging level
@@ -162,8 +176,11 @@ class AliCalorimeterUtils : public TObject {
   Bool_t             fRecalibration;         //  Switch on or off the recalibration
   TObjArray        * fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
   TObjArray        * fPHOSRecalibrationFactors;  // Array of histograms with map of recalibration factors, PHOS
-
-  ClassDef(AliCalorimeterUtils,2)
+  AliEMCALRecoUtils* fEMCALRecoUtils;        //  EMCAL utils for cluster rereconstruction
+  Bool_t             fRecalculatePosition;   // Recalculate cluster position
+  Bool_t             fCorrectELinearity  ;   // Correct cluster energy linearity
+  
+  ClassDef(AliCalorimeterUtils,3)
 } ;