1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // This analysis provides a new list of clusters to be used in other analysis
19 // Author: Gustavo Conesa Balbastre,
20 // Adapted from analysis class from Deepa Thomas
23 //_________________________________________________________________________
27 #include <TRefArray.h>
28 #include <TClonesArray.h>
30 #include <TGeoManager.h>
32 #include <TInterpreter.h>
35 // --- AliRoot Analysis Steering
36 #include "AliAnalysisTask.h"
37 #include "AliAnalysisManager.h"
38 #include "AliESDEvent.h"
39 #include "AliGeomManager.h"
40 #include "AliVCaloCells.h"
41 #include "AliAODCaloCluster.h"
42 #include "AliCDBManager.h"
43 #include "AliCDBStorage.h"
44 #include "AliCDBEntry.h"
46 #include "AliVEventHandler.h"
47 #include "AliAODInputHandler.h"
48 #include "AliOADBContainer.h"
49 #include "AliAODMCParticle.h"
52 #include "AliEMCALAfterBurnerUF.h"
53 #include "AliEMCALGeometry.h"
54 #include "AliEMCALClusterizerNxN.h"
55 #include "AliEMCALClusterizerv1.h"
56 #include "AliEMCALClusterizerv2.h"
57 #include "AliEMCALRecPoint.h"
58 #include "AliEMCALDigit.h"
59 #include "AliCaloCalibPedestal.h"
60 #include "AliEMCALCalibData.h"
62 #include "AliAnalysisTaskEMCALClusterize.h"
64 ClassImp(AliAnalysisTaskEMCALClusterize)
66 //______________________________________________________________________________
67 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
68 : AliAnalysisTaskSE(name)
70 , fGeom(0), fGeomName("")
71 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
72 , fCalibData(0), fPedestalData(0)
73 , fOCDBpath(""), fAccessOCDB(kFALSE)
74 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
75 , fRecParam(0), fClusterizer(0)
76 , fUnfolder(0), fJustUnfold(kFALSE)
77 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
78 , fFillAODFile(kFALSE), fFillAODHeader(0)
79 , fFillAODCaloCells(0), fRun(-1)
80 , fRecoUtils(0), fConfigName("")
82 , fCellLabels(), fCellSecondLabels(), fCellTime()
83 , fCellMatchdEta(), fCellMatchdPhi()
84 , fRecalibrateWithClusterTime(0)
85 , fMaxEvent(0), fDoTrackMatching(kFALSE)
86 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
87 , fRejectBelowThreshold(kFALSE)
88 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
89 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
90 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
91 , fCentralityClass(""), fSelectEMCALEvent(0)
92 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
93 , fSetCellMCLabelFromCluster(0)
94 , fRemapMCLabelForAODs(0)
99 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
103 fCentralityBin[0] = fCentralityBin[1]=-1;
107 //______________________________________________________________
108 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
109 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
111 , fGeom(0), fGeomName("")
112 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
113 , fCalibData(0), fPedestalData(0)
114 , fOCDBpath(""), fAccessOCDB(kFALSE)
115 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
116 , fRecParam(0), fClusterizer(0)
117 , fUnfolder(0), fJustUnfold(kFALSE)
118 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
119 , fFillAODFile(kFALSE), fFillAODHeader(0)
120 , fFillAODCaloCells(0), fRun(-1)
121 , fRecoUtils(0), fConfigName("")
122 , fOrgClusterCellId()
123 , fCellLabels(), fCellSecondLabels(), fCellTime()
124 , fCellMatchdEta(), fCellMatchdPhi()
125 , fRecalibrateWithClusterTime(0)
126 , fMaxEvent(0), fDoTrackMatching(kFALSE)
127 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
128 , fRejectBelowThreshold(kFALSE)
129 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
130 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
131 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
132 , fCentralityClass(""), fSelectEMCALEvent(0)
133 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
134 , fSetCellMCLabelFromCluster(0)
135 , fRemapMCLabelForAODs(0)
136 , fInputFromFilter(0)
140 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
144 fCentralityBin[0] = fCentralityBin[1]=-1;
149 //_______________________________________________________________
150 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
156 fDigitsArr->Clear("C");
162 fClusterArr->Delete();
168 fCaloClusterArr->Delete();
169 delete fCaloClusterArr;
172 if(fClusterizer) delete fClusterizer;
173 if(fUnfolder) delete fUnfolder;
174 if(fRecoUtils) delete fRecoUtils;
178 //_______________________________________________________
179 Bool_t AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL()
181 // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic
183 if(!fSelectEMCALEvent) return kTRUE; // accept
185 if(fEMCALEnergyCut <= 0) return kTRUE; // accept
187 Int_t nCluster = fEvent -> GetNumberOfCaloClusters();
188 AliVCaloCells * caloCell = fEvent -> GetEMCALCells();
189 Int_t bc = fEvent -> GetBunchCrossNumber();
191 for(Int_t icalo = 0; icalo < nCluster; icalo++)
193 AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
195 if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
196 fRecoUtils->IsGoodCluster(clus,fGeom,caloCell,bc))
200 printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Accept : E %2.2f > %2.2f, nCells %d > %d \n",
201 clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut);
209 printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Reject \n");
215 //_______________________________________________
216 void AliAnalysisTaskEMCALClusterize::AccessOADB()
218 // Set the AODB calibration, bad channels etc. parameters at least once
219 // alignment matrices from OADB done in SetGeometryMatrices
222 if(fOADBSet) return ;
224 Int_t runnumber = InputEvent()->GetRunNumber() ;
225 TString pass = GetPass();
227 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data());
229 Int_t nSM = fGeom->GetNumberOfSuperModules();
232 if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
234 AliOADBContainer *contBC=new AliOADBContainer("");
235 contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
237 TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
241 fRecoUtils->SwitchOnDistToBadChannelRecalculation();
242 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n");
244 for (Int_t i=0; i<nSM; ++i)
246 TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
251 hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
255 AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
259 hbm->SetDirectory(0);
260 fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
263 } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
266 // Energy Recalibration
267 if(fRecoUtils->IsRecalibrationOn())
269 AliOADBContainer *contRF=new AliOADBContainer("");
271 contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
273 TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
277 TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
281 TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
285 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n");
286 for (Int_t i=0; i<nSM; ++i)
288 TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
293 h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
297 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
303 fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
305 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array \n"); // array ok
306 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
307 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n"); // run number array ok
309 } // Recalibration on
311 // Energy Recalibration, apply on top of previous calibration factors
312 if(fRecoUtils->IsRunDepRecalibrationOn())
314 AliOADBContainer *contRFTD=new AliOADBContainer("");
316 contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
318 TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
320 //If it did not exist for this run, get closes one
323 AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
324 // let's get the closest runnumber instead then..
327 Int_t maxEntry = contRFTD->GetNumberOfEntries();
329 while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
334 Int_t closest = lower;
335 if ( (ic<maxEntry) &&
336 (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
340 AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
341 htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
346 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
348 for (Int_t ism=0; ism<nSM; ++ism)
350 for (Int_t icol=0; icol<48; ++icol)
352 for (Int_t irow=0; irow<24; ++irow)
354 Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
356 Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
357 factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
358 //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
359 // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
360 fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
364 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
365 } // Run by Run T calibration
367 // Time Recalibration
368 if(fRecoUtils->IsTimeRecalibrationOn())
370 AliOADBContainer *contTRF=new AliOADBContainer("");
372 contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
374 TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
378 TString passM = pass;
379 if(pass=="spc_calo") passM = "pass1";
380 TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
384 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Time Recalibrate EMCAL \n");
385 for (Int_t ibc = 0; ibc < 4; ++ibc)
387 TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
392 h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
396 AliError(Form("Could not load hAllTimeAvBC%d",ibc));
402 fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
403 } // bunch crossing loop
404 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
405 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n"); // run number array ok
407 } // Time recalibration on
409 // Parameters already set once, so do not it again
414 //_________________________________________________
415 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
417 //Access to OCDB stuff
419 fEvent = InputEvent();
422 Warning("AccessOCDB","Event not available!!!");
426 if (fEvent->GetRunNumber()==fRun)
428 fRun = fEvent->GetRunNumber();
430 if(DebugLevel() > 1 )
431 printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Begin");
433 AliCDBManager *cdb = AliCDBManager::Instance();
436 if (fOCDBpath.Length())
438 cdb->SetDefaultStorage(fOCDBpath.Data());
439 printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
442 cdb->SetRun(fEvent->GetRunNumber());
445 // EMCAL from RAW OCDB
446 if (fOCDBpath.Contains("alien:"))
448 cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
449 cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
452 TString path = cdb->GetDefaultStorage()->GetBaseFolder();
456 //Get calibration parameters
459 AliCDBEntry *entry = (AliCDBEntry*)
460 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
462 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
466 AliFatal("Calibration parameters not found in CDB!");
468 //Get calibration parameters
471 AliCDBEntry *entry = (AliCDBEntry*)
472 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
474 if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
478 AliFatal("Dead map not found in CDB!");
483 //_____________________________________________________
484 void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
486 // Get the input event, it can depend in embedded events what you want to get
487 // Also check if the quality of the event is good if not reject it
491 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
492 Int_t eventN = Entry();
493 if(aodIH) eventN = aodIH->GetReadEntry();
495 if (eventN > fMaxEvent)
498 //printf("Clusterizer --- Event %d-- \n",eventN);
500 //Check if input event are embedded events
501 //If so, take output event
502 if (aodIH && aodIH->GetMergeEvents())
506 if(!aodIH->GetMergeEMCALCells())
507 AliFatal("Events merged but not EMCAL cells, check analysis settings!");
511 printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
513 printf("\t InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
514 InputEvent()->GetEMCALCells()->GetNumberOfCells());
516 printf("\t MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
517 aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
519 for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
521 AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
522 if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E %f\n",icl,sigCluster->E());
525 printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
526 AODEvent()->GetEMCALCells()->GetNumberOfCells());
529 else if(fInputFromFilter)
531 //printf("Get Input From Filtered AOD\n");
536 fEvent = InputEvent();
537 if(fFillAODCaloCells) FillAODCaloCells();
538 if(fFillAODHeader) FillAODHeader();
543 Error("UserExec","Event not available");
547 //Process events if there is a high energy cluster
548 if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
550 //-------------------------------------------------------------------------------------
551 // Reject events if LED was firing, use only for LHC11a data
552 // Reject event if triggered by exotic cell and remove exotic cells if not triggered
553 //-------------------------------------------------------------------------------------
555 if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
557 if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
559 //-------------------------------------------------------------------------------------
560 // Set the cluster array in the event (output or input)
561 //-------------------------------------------------------------------------------------
565 //Magic line to write events to AOD filem put after event rejection
566 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
568 else if( !fOutputAODBranchSet )
570 // Create array and put it in the input event, if output AOD not selected, only once
571 InputEvent()->AddObject(fOutputAODBranch);
572 fOutputAODBranchSet = kTRUE;
573 printf("AliAnalysisTaskEMCALClusterize::UserExec() - Add AOD branch <%s> to input event\n",fOutputAODBranchName.Data());
578 //____________________________________________________
579 void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
581 // Recluster calocells, transform them into digits,
582 // feed the clusterizer with them and get new list of clusters
584 //In case of MC, first loop on the clusters and fill MC label to array
585 Int_t nClusters = fEvent->GetNumberOfCaloClusters();
586 Int_t nClustersOrg = 0;
588 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
589 if(aodIH && aodIH->GetEventToMerge()) //Embedding
590 nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
594 for (Int_t i = 0; i < nClusters; i++)
596 AliVCluster *clus = 0;
597 if(aodIH && aodIH->GetEventToMerge()) //Embedding
598 clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
600 clus = fEvent->GetCaloCluster(i);
606 Int_t label = clus->GetLabel();
608 //printf("Org cluster E %f, Time %e, Index = %d, ID %d, MC label %d\n", clus->E(), clus->GetTOF(),i, clus->GetID(),label );
609 //printf("Original list of labels from old cluster : \n");
610 //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
612 if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
613 UShort_t * index = clus->GetCellsAbsId() ;
614 for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
616 //printf("\t cell %d, MC label %d\n",index[icell],fEvent->GetEMCALCells()->GetCellMCLabel(index[icell]));
617 fOrgClusterCellId[index[icell]] = i;
618 fCellLabels[index[icell]] = label;
619 fCellSecondLabels[index[icell]] = label2;
620 fCellTime[index[icell]] = clus->GetTOF();
621 fCellMatchdEta[index[icell]] = clus->GetTrackDz();
622 fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
629 // Transform CaloCells into Digits
636 AliVCaloCells *cells = fEvent->GetEMCALCells();
638 Int_t bc = InputEvent()->GetBunchCrossNumber();
640 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
642 // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
643 id = cells->GetCellNumber(icell);
644 Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
646 // Do not include cells with too low energy, nor exotic cell
647 if( amp < fRecParam->GetMinECut() ||
648 time > fRecParam->GetTimeMax() ||
649 time < fRecParam->GetTimeMin() ) accept = kFALSE;
651 // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
652 if (fRecalibrateWithClusterTime)
654 time = fCellTime[id];
655 //printf("cell %d time org %f - ",id, time*1.e9);
656 fRecoUtils->RecalibrateCellTime(id,bc,time);
657 //printf("recal %f\n",time*1.e9);
661 if (accept && fRecoUtils->IsExoticCell(id,cells,bc))
666 if( DebugLevel() > 2 )
667 printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
668 id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
672 Int_t mcLabel = cells->GetMCLabel(icell);
673 //printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - cell %d, mc label %d\n",id,mcLabel);
675 //if(fCellLabels[id]!=mcLabel)printf("mcLabel %d - %d\n",mcLabel,fCellLabels[id]);
676 if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
677 else if( fSetCellMCLabelFromCluster == 0 && fRemapMCLabelForAODs) RemapMCLabelForAODs(mcLabel);
678 else mcLabel = -1; // found later
680 //printf("\t new label %d\n",mcLabel);
682 // Create the digit, put a fake primary deposited energy to trick the clusterizer
683 // when checking the most likely primary
685 Float_t efrac = cells->GetEFraction(icell);
687 //When checking the MC of digits, give weight to cells with embedded signal
688 if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
690 //printf("******* Cell %d, id %d, e %f, fraction %f, MC label %d, used MC label %d\n",icell,id,amp,cells->GetEFraction(icell),cells->GetMCLabel(icell),mcLabel);
692 new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, amp*efrac);
693 // Last parameter should be MC deposited energy, since it is not available, add just the cell amplitude so that
694 // we give more weight to the MC label of the cell with highest energy in the cluster
701 //-------------------------------------------------------------------------------------
702 //Do the clusterization
703 //-------------------------------------------------------------------------------------
705 fClusterizer->Digits2Clusters("");
707 //-------------------------------------------------------------------------------------
708 //Transform the recpoints into AliVClusters
709 //-------------------------------------------------------------------------------------
711 RecPoints2Clusters();
715 printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
719 if( DebugLevel() > 0 )
721 printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
723 if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
725 printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
726 fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
727 fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
732 //_____________________________________________________
733 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
735 // Take the event clusters and unfold them
737 AliVCaloCells *cells = fEvent->GetEMCALCells();
738 Double_t cellAmplitude = 0;
739 Double_t cellTime = 0;
740 Short_t cellNumber = 0;
741 Int_t cellMCLabel = 0;
742 Double_t cellEFrac = 0;
743 Int_t nClustersOrg = 0;
745 // Fill the array with the EMCAL clusters, copy them
746 for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
748 AliVCluster *clus = fEvent->GetCaloCluster(i);
751 //recalibrate/remove bad channels/etc if requested
752 if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
757 if(fRecoUtils->IsRecalibrationOn())
760 fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
763 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
765 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
768 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
769 fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
770 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
772 //Do not include bad channels found in analysis?
773 if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
774 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
777 cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
782 //Cast to ESD or AOD, needed to create the cluster array
783 AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
784 AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
788 fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
792 fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
795 Warning("UserExec()"," - Wrong CaloCluster type?");
802 fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
809 //_____________________________________________________
810 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells()
812 // Put calo cells in standard branch
813 AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
814 Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
816 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
817 aodEMcells.CreateContainer(nEMcell);
818 aodEMcells.SetType(AliVCaloCells::kEMCALCell);
819 Double_t calibFactor = 1.;
820 for (Int_t iCell = 0; iCell < nEMcell; iCell++)
822 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
823 fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
824 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
826 if(fRecoUtils->IsRecalibrationOn())
828 calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
831 if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
832 { //Channel is not declared as bad
833 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
834 eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
838 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
845 //__________________________________________________
846 void AliAnalysisTaskEMCALClusterize::FillAODHeader()
848 //Put event header information in standard AOD branch
850 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
851 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
855 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
857 AliAODHeader* header = AODEvent()->GetHeader();
858 header->SetRunNumber(fEvent->GetRunNumber());
862 TTree* tree = fInputHandler->GetTree();
865 TFile* file = tree->GetCurrentFile();
866 if (file) header->SetESDFileName(file->GetName());
869 else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
871 header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
872 header->SetOrbitNumber(fEvent->GetOrbitNumber());
873 header->SetPeriodNumber(fEvent->GetPeriodNumber());
874 header->SetEventType(fEvent->GetEventType());
877 if(fEvent->GetCentrality())
878 header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
880 header->SetCentrality(0);
883 header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
885 header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
887 header->SetTriggerMask(fEvent->GetTriggerMask());
888 header->SetTriggerCluster(fEvent->GetTriggerCluster());
892 header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
893 header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
894 header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
898 header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
899 header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
900 header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
903 header->SetMagneticField(fEvent->GetMagneticField());
904 //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
906 header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
907 header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
908 header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
909 header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
910 header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
912 Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
914 fEvent->GetDiamondCovXY(diamcov);
915 header->SetDiamond(diamxy,diamcov);
916 if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
917 else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
920 Int_t nVertices = 1 ;/* = prim. vtx*/;
921 Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
923 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
925 // Access to the AOD container of vertices
926 TClonesArray &vertices = *(AODEvent()->GetVertices());
929 // Add primary vertex. The primary tracks will be defined
930 // after the loops on the composite objects (V0, cascades, kinks)
931 fEvent->GetPrimaryVertex()->GetXYZ(pos);
935 esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
936 chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
940 aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
941 chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
944 AliAODVertex * primary = new(vertices[jVertices++])
945 AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
946 primary->SetName(fEvent->GetPrimaryVertex()->GetName());
947 primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
951 //___________________________________________________________
952 void AliAnalysisTaskEMCALClusterize::FillCaloClusterInEvent()
954 // Get the CaloClusters array, do some final calculations
955 // and put the clusters in the output or input event
956 // as a separate branch
958 //First recalculate track-matching for the new clusters
961 fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
963 //Put the new clusters in the AOD list
965 Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
967 for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
969 AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
971 newCluster->SetID(i);
973 // Correct cluster energy non linearity
974 newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
979 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
982 newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
984 printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
987 Float_t dR = 999., dZ = 999.;
988 fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
989 newCluster->SetTrackDistance(dR,dZ);
993 {// Assign previously assigned matched track in reco, very very rough
994 Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
995 newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
998 //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
999 //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1002 // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1003 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1005 new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
1007 if(DebugLevel() > 1 )
1008 printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f, mc label %d \n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel());
1012 fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1017 //_______________________________________________
1018 TString AliAnalysisTaskEMCALClusterize::GetPass()
1020 // Get passx from filename.
1022 if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1024 AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n");
1028 if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1030 AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n");
1034 TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1035 if (pass.Contains("ass1")) return TString("pass1");
1036 else if (pass.Contains("ass2")) return TString("pass2");
1037 else if (pass.Contains("ass3")) return TString("pass3");
1038 else if (pass.Contains("ass4")) return TString("pass4");
1039 else if (pass.Contains("ass5")) return TString("pass5");
1040 else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1041 else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1043 printf("AliAnalysisTaskEMCALClusterize::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n");
1044 return TString("pass1");
1047 // No condition fullfilled, give a default value
1048 printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
1053 //_________________________________________
1054 void AliAnalysisTaskEMCALClusterize::Init()
1056 //Init analysis with configuration macro if available
1059 if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
1061 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1063 if(!fRecParam) fRecParam = new AliEMCALRecParam;
1064 if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1066 if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1067 if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1068 if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1069 fRejectBelowThreshold = kFALSE;
1072 if(fCentralityClass == "") fCentralityClass = "V0M";
1074 if (fOCDBpath == "") fOCDBpath = "raw://" ;
1075 if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1077 if(gROOT->LoadMacro(fConfigName) >=0)
1079 printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
1080 AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1081 fGeomName = clus->fGeomName;
1082 fLoadGeomMatrices = clus->fLoadGeomMatrices;
1083 fOCDBpath = clus->fOCDBpath;
1084 fAccessOCDB = clus->fAccessOCDB;
1085 fRecParam = clus->fRecParam;
1086 fJustUnfold = clus->fJustUnfold;
1087 fFillAODFile = clus->fFillAODFile;
1088 fRecoUtils = clus->fRecoUtils;
1089 fConfigName = clus->fConfigName;
1090 fMaxEvent = clus->fMaxEvent;
1091 fDoTrackMatching = clus->fDoTrackMatching;
1092 fOutputAODBranchName = clus->fOutputAODBranchName;
1093 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1094 fCentralityClass = clus->fCentralityClass;
1095 fCentralityBin[0] = clus->fCentralityBin[0];
1096 fCentralityBin[1] = clus->fCentralityBin[1];
1101 //_______________________________________________________
1102 void AliAnalysisTaskEMCALClusterize::InitClusterization()
1104 //Select clusterization/unfolding algorithm and set all the needed parameters
1108 // init the unfolding afterburner
1110 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1114 //First init the clusterizer
1115 delete fClusterizer;
1116 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1117 fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
1118 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1119 fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
1120 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1122 fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
1123 fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1124 fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1128 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1131 //Now set the parameters
1132 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1133 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1134 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1135 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1136 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1137 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1138 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1139 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1140 fClusterizer->SetInputCalibrated ( kTRUE );
1141 fClusterizer->SetJustClusters ( kTRUE );
1143 // Initialize the cluster rec points and digits arrays and get them.
1144 fClusterizer->SetOutput(0);
1145 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1146 fDigitsArr = fClusterizer->GetDigits();
1148 //In case of unfolding after clusterization is requested, set the corresponding parameters
1149 if(fRecParam->GetUnfold())
1152 for (i = 0; i < 8; i++)
1154 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1155 }//end of loop over parameters
1157 for (i = 0; i < 3; i++)
1159 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1160 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1161 }//end of loop over parameters
1162 fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1163 fClusterizer->InitClusterUnfolding();
1168 //_________________________________________________
1169 void AliAnalysisTaskEMCALClusterize::InitGeometry()
1171 // Init geometry and set the geometry matrix, for the first event, skip the rest
1172 // Also set once the run dependent calibrations
1174 if(fGeomMatrixSet) return;
1176 Int_t runnumber = InputEvent()->GetRunNumber() ;
1181 if (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
1182 else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
1183 else fGeomName = "EMCAL_COMPLETE12SMV1";
1184 printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",
1185 fGeomName.Data(),runnumber);
1188 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1190 // Init geometry, I do not like much to do it like this ...
1191 if(fImportGeometryFromFile && !gGeoManager)
1193 if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1195 // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1196 if (runnumber < 140000 &&
1197 runnumber >= 100000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2010.root";
1198 if (runnumber >= 140000 &&
1199 runnumber < 171000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root";
1200 else fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1202 printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Import %s\n",fImportGeometryFilePath.Data());
1203 TGeoManager::Import(fImportGeometryFilePath) ;
1208 printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
1209 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1212 } // geometry pointer did not exist before
1214 if(fLoadGeomMatrices)
1216 printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n");
1218 // OADB if available
1219 AliOADBContainer emcGeoMat("AliEMCALgeo");
1220 emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1221 TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1224 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1227 if (!fGeomMatrix[mod]) // Get it from OADB
1230 printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n",
1231 mod,((TGeoHMatrix*) matEMCAL->At(mod)));
1232 //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1234 fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1237 if(fGeomMatrix[mod])
1239 if(DebugLevel() > 1)
1240 fGeomMatrix[mod]->Print();
1242 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1245 fGeomMatrixSet=kTRUE;
1249 else if(!gGeoManager)
1251 printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1252 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1253 if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1255 if(DebugLevel() > 1)
1256 Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
1260 if(DebugLevel() > 1)
1261 printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
1263 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1267 Error("InitGeometry"," - This event does not contain ESDs?");
1268 if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1272 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1274 if(DebugLevel() > 1)
1275 esd->GetEMCALMatrix(mod)->Print();
1277 if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1281 fGeomMatrixSet=kTRUE;
1284 }//Load matrices from Data
1288 //____________________________________________________
1289 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1292 // Check if event is exotic, get an exotic cell and compare with triggered patch
1293 // If there is a match remove event ... to be completed, filled with something provisional
1295 if(!fRemoveExoticEvents) return kFALSE;
1298 AliVCaloCells * cells = fEvent->GetEMCALCells();
1299 Float_t totCellE = 0;
1300 Int_t bc = InputEvent()->GetBunchCrossNumber();
1301 for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1305 Double_t tcell = 0 ;
1307 Int_t absID = cells->GetCellNumber(icell);
1308 Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1309 if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1312 // TString triggerclasses = event->GetFiredTriggerClasses();
1313 // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1314 // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1318 //printf("TotE cell %f\n",totCellE);
1319 if(totCellE < 1) return kTRUE;
1325 //________________________________________________________________
1326 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
1328 //Check if event is LED
1330 if(!fRemoveLEDEvents) return kFALSE;
1332 //check events of LHC11a period
1333 if(run < 146858 || run > 146860) return kFALSE ;
1335 // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1336 Int_t ncellsSM3 = 0;
1337 AliVCaloCells * cells = fEvent->GetEMCALCells();
1338 for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1340 if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1343 TString triggerclasses = fEvent->GetFiredTriggerClasses();
1345 Int_t ncellcut = 21;
1346 if(triggerclasses.Contains("EMC")) ncellcut = 35;
1348 if( ncellsSM3 >= ncellcut)
1350 printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d\n",(Int_t)Entry(),ncellsSM3);
1351 if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1359 //_______________________________________________________
1360 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters()
1362 // Restore clusters from recPoints
1363 // Cluster energy, global position, cells and their amplitude
1364 // fractions are restored
1367 for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1369 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1371 const Int_t ncells = recPoint->GetMultiplicity();
1372 Int_t ncellsTrue = 0;
1374 if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1376 // cells and their amplitude fractions
1377 UShort_t absIds[ncells];
1378 Double32_t ratios[ncells];
1380 //For later check embedding
1381 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1383 Float_t clusterE = 0;
1384 for (Int_t c = 0; c < ncells; c++)
1386 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1388 absIds[ncellsTrue] = digit->GetId();
1389 ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1391 // In case of unfolding, remove digits with unfolded energy too low
1394 if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1396 if(DebugLevel() > 1)
1398 printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
1399 recPoint->GetEnergiesList()[c],digit->GetAmplitude());
1407 //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1408 clusterE +=recPoint->GetEnergiesList()[c];
1410 // In case of embedding, fill ratio with amount of signal,
1411 if (aodIH && aodIH->GetMergeEvents())
1413 //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1414 AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1415 AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1417 Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1418 //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1419 Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1420 //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1422 if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1423 //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1429 }// cluster cell loop
1433 if (DebugLevel() > 1)
1434 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
1435 recPoint->GetEnergy(), ncells);
1439 //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1441 if(clusterE < fRecParam->GetClusteringThreshold())
1444 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
1451 // calculate new cluster position
1453 recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1454 recPoint->GetGlobalPosition(gpos);
1457 // create a new cluster
1459 (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1460 AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1462 clus->SetType(AliVCluster::kEMCALClusterv1);
1463 clus->SetE(clusterE);
1464 clus->SetPosition(g);
1465 clus->SetNCells(ncellsTrue);
1466 clus->SetCellsAbsId(absIds);
1467 clus->SetCellsAmplitudeFraction(ratios);
1468 clus->SetChi2(-1); //not yet implemented
1469 clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1470 clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1471 clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1473 if(ncells == ncellsTrue)
1475 Float_t elipAxis[2];
1476 recPoint->GetElipsAxis(elipAxis);
1477 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1478 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1479 clus->SetDispersion(recPoint->GetDispersion());
1481 else if(fSelectCell)
1483 // In case some cells rejected, in unfolding case, recalculate
1484 // shower shape parameters and position
1485 if(DebugLevel() > 1)
1486 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue);
1488 AliVCaloCells* cells = 0x0;
1489 if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1490 else cells = InputEvent()->GetEMCALCells();
1492 fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1493 fRecoUtils->RecalculateClusterPID(clus);
1494 fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1500 if ( fSetCellMCLabelFromCluster == 1 ) SetClustersMCLabelFrom2SelectedLabels(recPoint,clus) ;
1501 else if( fSetCellMCLabelFromCluster == 2 ) SetClustersMCLabelFromOriginalClusters(clus) ;
1504 // Normal case, trust what the clusterizer has found
1505 Int_t parentMult = 0;
1506 Int_t *parentList = recPoint->GetParents(parentMult);
1507 clus->SetLabel(parentList, parentMult);
1508 // printf("Label list : ");
1509 // for(Int_t ilabel = 0; ilabel < parentMult; ilabel++ ) printf(" %d ",parentList[ilabel]);
1517 //___________________________________________________________________________
1518 void AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs(Int_t & label)
1520 // MC label for Cells not remapped after ESD filtering, do it here.
1522 if(label < 0) return ;
1524 AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
1527 TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1530 if(label < arr->GetEntriesFast())
1532 AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1533 if(!particle) return ;
1535 if(label == particle->Label()) return ; // label already OK
1536 //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
1538 //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
1540 // loop on the particles list and check if there is one with the same label
1541 for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
1543 AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1544 if(!particle) continue ;
1546 if(label == particle->Label())
1549 //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
1556 //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
1560 //________________________________________________
1561 void AliAnalysisTaskEMCALClusterize::ResetArrays()
1563 // Reset arrays containing information for all possible cells
1564 for(Int_t j = 0; j < 12672; j++)
1566 fOrgClusterCellId[j] = -1;
1567 fCellLabels[j] = -1;
1568 fCellSecondLabels[j] = -1;
1570 fCellMatchdEta[j] = -999;
1571 fCellMatchdPhi[j] = -999;
1575 //_____________________________________________________________________________________________________
1576 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint * recPoint,
1577 AliAODCaloCluster * clus)
1579 // Set the cluster MC label, the digizer was filled with most likely MC label for all cells in original cluster
1580 // Now check the second most likely MC label and add it to the new cluster
1582 Int_t parentMult = 0;
1583 Int_t *parentList = recPoint->GetParents(parentMult);
1584 clus->SetLabel(parentList, parentMult);
1586 //Write the second major contributor to each MC cluster.
1588 for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1591 Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1594 iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
1595 for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1597 if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
1598 if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
1602 Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1603 for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1605 newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1608 newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1609 clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1610 delete [] newLabelArray;
1612 }//positive cell number
1616 //___________________________________________________________________________________________________
1617 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster * clus)
1619 // Get the original clusters that contribute to the new cluster, assign the labels of such clusters
1620 // to the new cluster.
1621 // Only approximatedly valid when output are V1 clusters, handle with care
1623 TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1626 Int_t nLabTotOrg = 0;
1630 AliVEvent * event = fEvent;
1632 //In case of embedding the MC signal is in the added event
1633 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1634 if(aodIH && aodIH->GetEventToMerge()) //Embedding
1635 event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
1638 //Find the clusters that originally had the cells
1639 for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1641 Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1645 Int_t idCluster = fOrgClusterCellId[idCell];
1648 for(Int_t icl =0; icl < nClu; icl++)
1650 if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
1651 if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
1652 // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
1653 // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
1655 if( set && idCluster >= 0)
1657 clArray.SetAt(idCluster,nClu++);
1658 //printf("******** idCluster %d \n",idCluster);
1659 nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
1661 //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
1663 //Search highest E cluster
1664 AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1665 //printf("\t E %f\n",clOrg->E());
1666 if(emax < clOrg->E())
1675 if(nClu==0 || nLabTotOrg == 0)
1677 //if(clus->E() > 0.25) printf("AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters() - Check: N org clusters %d, n tot labels %d, cluster E %f, n cells %d\n",nClu,nLabTotOrg,clus->E(), clus->GetNCells());
1678 //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
1682 // Put the first in the list the cluster with highest energy
1683 if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1685 Int_t maxIndex = -1;
1686 Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1687 for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1689 if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1692 if(firstCluster >=0 && idMax >=0)
1694 clArray.SetAt(idMax,0);
1695 clArray.SetAt(firstCluster,maxIndex);
1699 // Get the labels list in the original clusters, assign all to the new cluster
1700 TArrayI clMCArray(nLabTotOrg) ;
1704 for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1706 Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
1707 //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
1708 AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1709 Int_t nLab = clOrg->GetNLabels();
1711 for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
1713 Int_t lab = clOrg->GetLabelAt(iLab) ;
1717 //printf("\t \t Set Label %d \n", lab);
1718 for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1720 if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
1721 //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
1722 // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
1724 if( set ) clMCArray.SetAt(lab,nLabTot++);
1729 // Set the final list of labels
1731 clus->SetLabel(clMCArray.GetArray(), nLabTot);
1733 // printf("Final list of labels for new cluster : \n");
1734 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
1736 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
1737 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
1738 // printf(" org %d ",label);
1739 // RemapMCLabelForAODs(label);
1740 // printf(" new %d \n",label);
1742 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
1746 //____________________________________________________________
1747 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1749 // Init geometry, create list of output clusters
1752 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1754 if(fOutputAODBranchName.Length()==0)
1756 fOutputAODBranchName = "newEMCALClustersArray";
1757 printf("Cluster branch name not set, set it to newEMCALClustersArray \n");
1760 fOutputAODBranch->SetName(fOutputAODBranchName);
1764 //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1767 //fOutputAODBranch->SetOwner(kFALSE);
1769 AddAODBranch("TClonesArray", &fOutputAODBranch);
1774 //_______________________________________________________
1775 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
1777 // Do clusterization event by event, execute different steps
1778 // 1) Do some checks on the kind of events (ESD, AOD) or if some filtering is needed, initializations
1779 // load things and clear arrays
1780 // 2) Clusterize a) just unfolding existing clusters (fJustUnfold)
1781 // b) recluster cells
1782 // + convert cells into digits (calibrating them)
1783 // + recluster digits into recPoints with the chosen clusterizer (V1, V1+Unfold,V2, NxN)
1784 // with methods in AliEMCALClusterizer
1785 // + transform recPoints into CaloClusters
1786 // 3) Do some final calculations in the found clusters (track-matching) and put them in an AOD branch
1792 //Remove the contents of AOD branch output list set in the previous event
1793 fOutputAODBranch->Clear("C");
1797 //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
1798 //If we need a centrality bin, we select only those events in the corresponding bin.
1799 if( GetCentrality() && fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1801 Float_t cen = GetEventCentrality();
1802 if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
1805 // intermediate array with new clusters : init the array only once or clear from previous event
1806 if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
1807 else fCaloClusterArr->Delete();//Clear("C"); it leaks?
1809 InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
1811 //Get the event, do some checks and settings
1812 CheckAndGetEvent() ;
1816 if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1820 //Init pointers, geometry, clusterizer, ocdb, aodb
1823 if(fAccessOCDB) AccessOCDB();
1824 if(fAccessOADB) AccessOADB(); // only once
1826 InitClusterization();
1832 if (fJustUnfold) ClusterUnfolding();
1833 else ClusterizeCells() ;
1838 FillCaloClusterInEvent();