EMCAL Tender Supply added.
[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliEMCALTenderSupply.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  EMCAL tender, apply corrections to EMCAl clusters                        //
20 //  and do track matching                                                    //                                                                           
21 //  Author : Deepa Thomas (Utrecht University)                               //                      
22 //                                                                           //
23 ///////////////////////////////////////////////////////////////////////////////
24
25 #include "TROOT.h"
26 #include "TFile.h"
27 #include "TAlienFile.h"
28 #include "TGrid.h"
29 #include "TTree.h"
30 #include "TInterpreter.h"
31 #include "TObjArray.h"
32
33 #include <AliLog.h>
34 #include <AliESDEvent.h>
35 #include <AliAnalysisManager.h>
36 #include <AliTender.h>
37 #include "AliOADBContainer.h"
38 #include "AliMagF.h"
39 #include "TGeoGlobalMagField.h"
40
41 #include "AliESDCaloCluster.h"
42 #include "AliEMCALTenderSupply.h"
43
44 //EMCAL
45 #include "AliEMCALGeometry.h"
46 #include "AliEMCALRecoUtils.h"
47
48 AliEMCALTenderSupply::AliEMCALTenderSupply() :
49   AliTenderSupply()
50  ,fEMCALGeo(0x0)
51  ,fEMCALGeoName("EMCAL_FIRSTYEARV1")
52  ,fEMCALRecoUtils(new AliEMCALRecoUtils)
53  ,fConfigName("")
54  ,fDebugLevel(0)
55  ,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection) 
56  ,fNonLinearThreshold(30)       
57  ,fReCalib(kFALSE)      
58  ,fRecalClusPos(kFALSE)
59  ,fFiducial(kFALSE) 
60  ,fNCellsFromEMCALBorder(1)     
61  ,fRecalDistToBadChannels(kFALSE)       
62  ,fInputTree(0) 
63  ,fInputFile(0)
64  ,fFilepass(0) 
65  ,fMass(0.139)
66  ,fStep(1)
67  ,fRcut(15)     
68 {
69   //
70   // default ctor
71   //
72         for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
73 }
74
75 //_____________________________________________________
76 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
77   AliTenderSupply(name,tender)
78  ,fEMCALGeo(0x0)
79  ,fEMCALGeoName("EMCAL_FIRSTYEARV1")
80  ,fEMCALRecoUtils(new AliEMCALRecoUtils)
81  ,fConfigName("") 
82  ,fDebugLevel(0)
83  ,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection)              
84  ,fNonLinearThreshold(30)       
85  ,fReCalib(kFALSE)      
86  ,fRecalClusPos(kFALSE)
87  ,fFiducial(kFALSE) 
88  ,fNCellsFromEMCALBorder(1)     
89  ,fRecalDistToBadChannels(kFALSE)       
90  ,fInputTree(0) 
91  ,fInputFile(0)
92  ,fFilepass(0)
93  ,fMass(0.139)
94  ,fStep(1)
95  ,fRcut(15)     
96 {
97   //
98   // named ctor
99   //
100         for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
101 }
102
103 //_____________________________________________________
104 AliEMCALTenderSupply::~AliEMCALTenderSupply()
105 {
106    //Destructor
107    delete fEMCALGeo;
108    delete fEMCALRecoUtils;
109    delete fInputTree;
110    delete fInputFile;
111 }
112
113 //_____________________________________________________
114 void AliEMCALTenderSupply::Init()
115 {
116   //
117   // Initialise EMCAL tender
118   //
119
120   if (fDebugLevel>0) AliInfo("Init EMCAL Tender supply\n");     
121
122   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
123  
124   fInputTree = mgr->GetTree();
125      
126   if(gROOT->LoadMacro(fConfigName) >=0){
127     AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()");
128     fDebugLevel         = tender->fDebugLevel;
129     fEMCALGeoName       = tender->fEMCALGeoName; 
130     fEMCALRecoUtils     = tender->fEMCALRecoUtils; 
131     fConfigName         = tender->fConfigName;
132     fNonLinearFunc      = tender->fNonLinearFunc;
133     fNonLinearThreshold = tender->fNonLinearThreshold;
134     fReCalib            = tender->fReCalib;
135     fRecalClusPos       = tender->fRecalClusPos;
136     fFiducial             = tender->fFiducial;
137     fNCellsFromEMCALBorder = tender->fNCellsFromEMCALBorder;
138     fRecalDistToBadChannels = tender->fRecalDistToBadChannels;    
139     fMass               = tender->fMass;
140     fStep               = tender->fStep;
141     fRcut               = tender->fRcut;   
142   }
143   
144   // Init goemetry      
145   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
146
147   fEMCALRecoUtils = new AliEMCALRecoUtils();
148
149   //Initialising Non linearity parameters
150   fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
151
152 }
153
154 //_____________________________________________________
155 void AliEMCALTenderSupply::ProcessEvent()
156 {
157   //Event loop  
158   AliESDEvent *event=fTender->GetEvent();
159   if (!event) return;
160
161   if(fTender->RunChanged()){ 
162         //Initialising parameters once per run number
163         if (!InitBadChannels()) return;
164         if (fReCalib){ if (!InitRecalibCluster()) return;}
165         if (fRecalClusPos){ if (!InitMisalignMatrix()) return;}
166   }
167
168   AliESDCaloCells *cells= event->GetEMCALCells();
169
170   //------------Good clusters------------
171  TClonesArray *clusArr;
172
173  clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
174  if(!clusArr) clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
175  if(!clusArr) return;
176
177  for (Int_t icluster=0; icluster < clusArr->GetEntries(); icluster++ ){
178          AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
179          if(!clust) continue;
180          if (!clust->IsEMCAL()) continue;
181          if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells() )){
182                  delete clusArr->RemoveAt(icluster);
183                  continue;
184          }
185          if(fFiducial){
186                  if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
187                         delete clusArr->RemoveAt(icluster);
188                         continue;
189                  }
190          }
191          fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
192          if(fRecalDistToBadChannels) fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);  
193          if(fReCalib) fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
194          if(fRecalClusPos) fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
195          //fEMCALRecoUtils->SetTimeDependentCorrections(event->GetRunNumber());
196  }
197  clusArr->Compress();
198
199  //--------Track Matching--------------
200
201  //set magnetic field
202  Double_t  magF = event->GetMagneticField();
203  Double_t magSign = 1.0;
204  if(magF<0)magSign = -1.0;
205   
206  if (!TGeoGlobalMagField::Instance()->GetField()) 
207           {
208                   AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG);
209                   TGeoGlobalMagField::Instance()->SetField(field);
210           }
211
212   //Setting mass, step size and residual cut  
213   fEMCALRecoUtils->SetCutR(fRcut);
214   fEMCALRecoUtils->SetMass(fMass);
215   fEMCALRecoUtils->SetStep(fStep);
216    
217   fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
218
219   Int_t NTracks = event->GetNumberOfTracks();
220
221   if(NTracks>0){
222         SetClusterMatchedToTrack(event);
223         SetTracksMatchedToCluster(event);
224   }
225 }
226
227 //_____________________________________________________
228 void AliEMCALTenderSupply::SetClusterMatchedToTrack(AliESDEvent *event)
229 {
230   //checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track. 
231
232   Int_t matchClusIndex = -1;
233   Int_t nTracks = event->GetNumberOfTracks();
234
235   // Track loop 
236   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) 
237   {
238            AliESDtrack* track = event->GetTrack(iTrack);
239            if (!track) {
240                    printf("ERROR: Could not receive track %d\n", iTrack);
241                    continue;
242            }
243
244            matchClusIndex = fEMCALRecoUtils->GetMatchedClusterIndex(iTrack);               
245            track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
246   }
247
248     if (fDebugLevel>1) AliInfo("Track matched to closest cluster\n");   
249 }
250
251 //_____________________________________________________
252 void AliEMCALTenderSupply::SetTracksMatchedToCluster(AliESDEvent *event)
253 {
254         //checks if EMC clusters are matched to ESD track.
255         //Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster
256         
257         Int_t matchTrackIndex = -1;
258         Int_t nTracks = event->GetNumberOfTracks();
259           
260         // Cluster loop
261         for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); iClus++)
262         {
263                 AliESDCaloCluster * cluster = event->GetCaloCluster(iClus);
264                 if (!cluster->IsEMCAL()) continue;
265                 
266                 Int_t nMatched =0;
267                 TArrayI arrayTrackMatched(nTracks);
268
269                 //get the closest track matched to the cluster
270                 matchTrackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(iClus); 
271                 arrayTrackMatched[nMatched] = matchTrackIndex;
272                 nMatched++;
273
274                 //get all other tracks matched to the cluster
275                 //track loop
276                 for(Int_t iTrk=0; iTrk<nTracks; iTrk++){
277                         AliESDtrack* track = event->GetTrack(iTrk);
278                         if(iTrk == matchTrackIndex) continue;
279
280                         if(track->GetEMCALcluster() == iClus){
281                                 arrayTrackMatched[nMatched] = iTrk;
282                                 nMatched++;
283                         }
284                 }
285                 arrayTrackMatched.Set(nMatched);
286                 cluster->AddTracksMatched(arrayTrackMatched);
287         
288         }
289     if (fDebugLevel>1) AliInfo("Cluster matched to tracks \n"); 
290 }
291
292 //_____________________________________________________
293 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
294 {
295   //Initialising Misalignment matrix
296         
297   AliESDEvent *event=fTender->GetEvent();
298   if (!event) return kFALSE;
299
300   Int_t runGM = 0;
301   runGM = event->GetRunNumber();
302         
303   if (fDebugLevel>0) AliInfo("Initialising Misalignment matrix \n");    
304
305   fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
306         
307   if(runGM <=140000){ //2010 data
308
309                 Double_t rotationMatrix[4][9] = {{-0.014587, -0.999892, -0.002031, 0.999892, -0.014591,  0.001979, -0.002009, -0.002002,  0.999996},
310                                                  {-0.014587,  0.999892,  0.002031, 0.999892,  0.014591, -0.001979, -0.002009,  0.002002, -0.999996},
311                                                  {-0.345864, -0.938278, -0.003412, 0.938276, -0.345874,  0.003010, -0.004004, -0.002161,  0.999990},
312                                                  {-0.345861,  0.938280,  0.003412, 0.938276,  0.345874, -0.003010, -0.004004,  0.002161, -0.999990}};
313   
314                 Double_t translationMatrix[4][3] = {{0.351659,    447.576446,  176.269742},
315                                                   {1.062577,    446.893974, -173.728870},
316                                                   {-154.213287, 419.306156,  176.753692},
317                                                   {-153.018950, 418.623681, -173.243605}};
318
319                 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
320                         {
321 //                              if(DebugLevel() > 1)  fEMCALMatrix[mod]->Print();
322                                 fEMCALMatrix[mod]->SetRotation(rotationMatrix[mod]);
323                                 fEMCALMatrix[mod]->SetTranslation(translationMatrix[mod]);              
324                                 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod); 
325                 }
326         }
327         
328
329   if(runGM>140000){ // 2011 data
330         AliOADBContainer emcalgeoCont(Form("emcal2011"));
331         emcalgeoCont.InitFromFile(Form("BetaGood.root"),Form("AliEMCALgeo"));
332
333         TObjArray *mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS");
334         
335         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
336                 fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
337                 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod); 
338
339                 fEMCALMatrix[mod]->Print();
340         }
341   }
342         
343  return kTRUE;
344 }
345
346 //_____________________________________________________
347 Bool_t AliEMCALTenderSupply::InitBadChannels()
348 {
349   //Initialising Bad channel maps
350         
351   AliESDEvent *event=fTender->GetEvent();
352   if (!event) return kFALSE;
353
354   Int_t fRunBC= 0;
355   fRunBC = event->GetRunNumber();
356
357   if (fDebugLevel>0) AliInfo("Initialising Bad channel map \n");        
358    
359   if(fFiducial){
360         fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
361         fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
362   }
363
364   fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
365   if(fRecalDistToBadChannels) fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
366
367   TFile *fbad;
368   //2010 
369   if(fRunBC <=140000){
370         fbad = new TFile("BadChannels.root","read");
371          if (fbad->IsZombie()){
372                 if (fDebugLevel>1) AliInfo("Connecting to alien to get BadChannels.root \n");   
373                 TGrid::Connect("alien://");
374                 fbad = TFile::Open("alien:///alice/cern.ch/user/g/gconesab/BadChannelsDB/BadChannels.root");
375         }
376         
377         TH2I * hb0 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod0");
378         TH2I * hb1 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod1");
379         TH2I * hb2 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod2");
380         TH2I * hb3 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod3"); 
381         fEMCALRecoUtils->SetEMCALChannelStatusMap(0,hb0);
382         fEMCALRecoUtils->SetEMCALChannelStatusMap(1,hb1);
383         fEMCALRecoUtils->SetEMCALChannelStatusMap(2,hb2);
384         fEMCALRecoUtils->SetEMCALChannelStatusMap(3,hb3); 
385   }
386
387   //2011
388   if(fRunBC>140000){
389             
390     const Int_t nTowers=89;
391     Int_t hotChannels[nTowers]={74, 103, 152, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 368, 369, 370, 371, 372, 373, 374,375, 376, 377, 378, 379, 380, 381, 382, 383, 917, 1275, 1288, 1519, 1595, 1860, 1967, 2022, 2026, 2047, 2117, 2298, 2540, 2776, 3135, 3764, 6095, 6111, 6481, 6592, 6800, 6801, 6802, 6803, 6804, 6805, 6806, 6807, 6808, 6809, 6810, 6811, 6812, 6813, 6814, 6815, 7371, 7425, 7430, 7457, 7491, 7709, 8352, 8353, 8356, 8357, 8808, 8810, 8812, 8814, 9056, 9769, 9815, 9837};
392
393     Int_t nSupMod=-1, nModule=-1, nIphi=-1, nIeta=-1, iphi=-1, ieta=-1;
394     for(Int_t i=0; i<nTowers; i++)
395             {
396                     fEMCALGeo->GetCellIndex(hotChannels[i],nSupMod,nModule,nIphi,nIeta);
397                     
398                     fEMCALGeo->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta,iphi,ieta);
399                     fEMCALRecoUtils->SetEMCALChannelStatus(nSupMod, ieta, iphi);
400             }
401
402    }
403
404    return kTRUE;
405 }
406
407 //_____________________________________________________
408 Bool_t AliEMCALTenderSupply::InitRecalibCluster()
409 {
410   //Initialising Recalibration Factors
411         
412   AliESDEvent *event=fTender->GetEvent();
413   if (!event) return kFALSE;
414
415   if(fInputTree){ 
416      fInputFile = fInputTree->GetCurrentFile();
417      if(fInputFile){
418             const char *fileName = fInputFile->GetName();
419             TString FileName = TString(fileName);
420             if     (FileName.Contains("pass1")) fFilepass = TString("pass1");
421             else if(FileName.Contains("pass2")) fFilepass = TString("pass2");
422             else if(FileName.Contains("pass3")) fFilepass = TString("pass3");
423             else AliError("pass number not found");
424      }
425       else AliError("File not found");
426   }
427   else AliError("Tree not found");
428 //        else {cout << "Tree not found " <<endl; return kFALSE;}
429
430 //      if (event->GetRunNumber()==runRC)
431 //          return kTRUE;
432
433   Int_t runRC= 0;           
434   runRC = event->GetRunNumber();
435
436   if (fDebugLevel>0) AliInfo("Initialising Recalibration factors \n");  
437
438    fEMCALRecoUtils->SwitchOnRecalibration();
439
440    TFile* fRecalib;
441    fRecalib = new TFile("RecalibrationFactors.root","read");
442    if (fRecalib->IsZombie()){
443         if (fDebugLevel>1) AliInfo("Connecting to alien to get RecalibrationFactors.root \n");  
444
445         if(     (runRC >= 114737 && runRC <= 117223 && (fFilepass = "pass2")) || //LHC10b pass2 
446                 (runRC >= 118503 && runRC <= 121040 && ((fFilepass = "pass2")||(fFilepass = "pass3"))) || //LHC10c pass2, LHC10c pass3
447                 (runRC >= 122195 && runRC <= 126437 && (fFilepass = "pass1")) || //LHC10d pass1
448                 (runRC >= 127712 && runRC <= 130850 && (fFilepass = "pass1")) || //LHC10e pass1
449                 (runRC >= 133004 && runRC < 134657  && (fFilepass = "pass1")))// LHC10f pass1 <134657
450                 {
451                         TGrid::Connect("alien://");
452                         fRecalib = TFile::Open("alien:///alice/cern.ch/user/g/gconesab/RecalDB/summer_december_2010/RecalibrationFactors.root");
453                 }
454
455                 else if((runRC >= 122195 && runRC <= 126437 && (fFilepass = "pass2")) || //LHC10d pass2
456                         (runRC >= 134657 && runRC <= 135029 && (fFilepass = "pass1")) || //LHC10f pass1 >= 134657
457                         (runRC >= 135654 && runRC <= 136377 && (fFilepass = "pass1")) || //LHC10g pass1
458                         (runRC >= 136851 && runRC < 137231  && (fFilepass = "pass1"))) //LHC10h pass1 untill christmas
459                 {
460                         TGrid::Connect("alien://");
461                         fRecalib = TFile::Open("alien:///alice/cern.ch/user/g/gconesab/RecalDB/december2010/RecalibrationFactors.root");
462                 }
463         
464                 else if(runRC >= 137231 && runRC <= 139517 && (fFilepass = "pass1")) //LHC10h pass1 from christmas
465                 {
466                         TGrid::Connect("alien://");
467                         fRecalib = TFile::Open("alien:///alice/cern.ch/user/g/gconesab/RecalDB/summer2010/RecalibrationFactors.root");
468                 }
469                 else {
470                         AliError("Run number not found OR Run number of 2011 for which recalibration not avaliable yet");
471                         return kFALSE;
472                 }
473     }
474
475     TH2F * r0 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM0");
476     TH2F * r1 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM1");
477     TH2F * r2 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM2");
478     TH2F * r3 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM3");
479     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(0,r0);
480     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(1,r1);
481     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(2,r2);
482     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(3,r3); 
483
484     return kTRUE;
485 }
486