X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerv1.cxx;h=8557f4754bdc459f157d140e356fb733b33b469e;hb=803d1ab0e2c9c9a8116c05954c3d3cdcf1f33759;hp=9c6bad43087cdbdaaad01d164b91bb690b1f7af9;hpb=11f9c5ffb14c571666639e9c11ac07d6d505a130;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 9c6bad43087..8557f4754bd 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -14,16 +14,7 @@ **************************************************************************/ /* $Id$ */ -/* $Log: - 1 October 2000. Yuri Kharlov: - AreNeighbours() - PPSD upper layer is considered if number of layers>1 - 18 October 2000. Yuri Kharlov: - AliEMCALClusterizerv1() - CPV clusterizing parameters added - MakeClusters() - After first PPSD digit remove EMC digits only once -*/ + //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction // of new IO (à la PHOS) @@ -68,15 +59,13 @@ // --- AliRoot header files --- - +#include "AliEMCALGetter.h" #include "AliEMCALClusterizerv1.h" +#include "AliEMCALTowerRecPoint.h" #include "AliEMCALDigit.h" #include "AliEMCALDigitizer.h" -#include "AliEMCALTowerRecPoint.h" #include "AliEMCAL.h" -#include "AliEMCALGetter.h" #include "AliEMCALGeometry.h" -#include "AliRun.h" ClassImp(AliEMCALClusterizerv1) @@ -90,8 +79,8 @@ ClassImp(AliEMCALClusterizerv1) } //____________________________________________________________________________ -AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit) -:AliEMCALClusterizer(headerFile, name, toSplit) +AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName) +:AliEMCALClusterizer(alirunFileName, eventFolderName) { // ctor with the indication of the file where header Tree and digits Tree are stored @@ -105,26 +94,30 @@ AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* AliEMCALClusterizerv1::~AliEMCALClusterizerv1() { // dtor - fSplitFile = 0 ; } //____________________________________________________________________________ const TString AliEMCALClusterizerv1::BranchName() const { - TString branchName(GetName() ) ; - branchName.Remove(branchName.Index(Version())-1) ; - return branchName ; + return GetName(); + } //____________________________________________________________________________ -Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Bool_t inpresho) const +Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t where) const { //To be replased later by the method, reading individual parameters from the database - if ( inpresho ) // calibrate as pre shower - return -fADCpedestalPreSho + amp * fADCchannelPreSho ; - else //calibrate as tower - return -fADCpedestalTower + amp * fADCchannelTower ; + // where = 0 == PRE ; where = 1 == ECAL ; where = 2 == HCAL + if ( where == 0 ) // calibrate as PRE section + return -fADCpedestalPRE + amp * fADCchannelPRE ; + else if (where == 1) //calibrate as ECA section + return -fADCpedestalECA + amp * fADCchannelECA ; + else if (where == 2) //calibrate as HCA section + return -fADCpedestalHCA + amp * fADCchannelHCA ; + else + Fatal("Calibrate", "Something went wrong!") ; + return -9999999. ; } //____________________________________________________________________________ @@ -132,18 +125,14 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) { // Steering method - if( strcmp(GetName(), "")== 0 ) - Init() ; - if(strstr(option,"tim")) gBenchmark->Start("EMCALClusterizer"); if(strstr(option,"print")) Print("") ; - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; - if(gime->BranchExists("RecPoints")) - return ; + AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + Int_t nevents = gime->MaxEvent() ; Int_t ievent ; @@ -151,13 +140,12 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) gime->Event(ievent,"D") ; - if(ievent == 0) - GetCalibrationParameters() ; + GetCalibrationParameters() ; - fNumberOfTowerClusters = fNumberOfPreShoClusters = 0 ; + fNumberOfPREClusters = fNumberOfECAClusters = fNumberOfHCAClusters = 0 ; MakeClusters() ; - + if(fToUnfold) MakeUnfolding() ; @@ -166,17 +154,19 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) if(strstr(option,"deb")) PrintRecPoints(option) ; - //increment the total number of digits per run - fRecPointsInRun += gime->TowerRecPoints()->GetEntriesFast() ; - fRecPointsInRun += gime->PreShowerRecPoints()->GetEntriesFast() ; - } + //increment the total number of recpoints per run + fRecPointsInRun += gime->PRERecPoints()->GetEntriesFast() ; + fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ; + fRecPointsInRun += gime->HCARecPoints()->GetEntriesFast() ; + } + Unload(); + if(strstr(option,"tim")){ gBenchmark->Stop("EMCALClusterizer"); Info("Exec", "took %f seconds for Clusterizing %f seconds per event", gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nevents ) ; } - } //____________________________________________________________________________ @@ -187,7 +177,7 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDig // The initial values for fitting procedure are set equal to the positions of local maxima. // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; + AliEMCALGetter * gime = AliEMCALGetter::Instance() ; TClonesArray * digits = gime->Digits() ; @@ -274,14 +264,20 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDig //____________________________________________________________________________ void AliEMCALClusterizerv1::GetCalibrationParameters() { - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; - const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ; + AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + + if ( !gime->Digitizer() ) + gime->LoadDigitizer(); + AliEMCALDigitizer * dig = gime->Digitizer(); + + fADCchannelPRE = dig->GetPREchannel() ; + fADCpedestalPRE = dig->GetPREpedestal() ; - fADCchannelTower = dig->GetTowerchannel() ; - fADCpedestalTower = dig->GetTowerpedestal(); + fADCchannelECA = dig->GetECAchannel() ; + fADCpedestalECA = dig->GetECApedestal(); - fADCchannelPreSho = dig->GetPreShochannel() ; - fADCpedestalPreSho = dig->GetPreShopedestal() ; + fADCchannelHCA = dig->GetHCAchannel() ; + fADCpedestalHCA = dig->GetHCApedestal(); } //____________________________________________________________________________ @@ -290,75 +286,39 @@ void AliEMCALClusterizerv1::Init() // Make all memory allocations which can not be done in default constructor. // Attach the Clusterizer task to the list of EMCAL tasks - if ( strcmp(GetTitle(), "") == 0 ) - SetTitle("galice.root") ; + AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data()); - TString branchname = GetName() ; - branchname.Remove(branchname.Index(Version())-1) ; - - AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ; - if ( gime == 0 ) { - Error("Init", "Could not obtain the Getter object !" ) ; - return ; - } - - fSplitFile = 0 ; - if(fToSplit){ - // construct the name of the file as /path/EMCAL.SDigits.root - //First - extract full path if necessary - TString fileName(GetTitle()) ; - Ssiz_t islash = fileName.Last('/') ; - if(islash(gROOT->GetFile(fileName.Data())); - if(!fSplitFile) - fSplitFile = TFile::Open(fileName.Data(),"update") ; - } + AliEMCALGeometry * geom = gime->EMCALGeometry() ; - const AliEMCALGeometry * geom = gime->EMCALGeometry() ; fNTowers = geom->GetNZ() * geom->GetNPhi() ; if(!gMinuit) gMinuit = new TMinuit(100) ; - gime->PostClusterizer(this) ; - gime->PostRecPoints(branchname ) ; - + if ( !gime->Clusterizer() ) + gime->PostClusterizer(this); } //____________________________________________________________________________ void AliEMCALClusterizerv1::InitParameters() { - fNumberOfPreShoClusters = fNumberOfTowerClusters = 0 ; - fPreShoClusteringThreshold = 0.0001; - fTowerClusteringThreshold = 0.2; - fTowerLocMaxCut = 0.03 ; - fPreShoLocMaxCut = 0.03 ; + fNumberOfPREClusters = fNumberOfECAClusters = fNumberOfHCAClusters = 0 ; + fPREClusteringThreshold = 0.0001; // must be adjusted according to the noise leve set by digitizer + fECAClusteringThreshold = 0.0045; // must be adjusted according to the noise leve set by digitizer + fHCAClusteringThreshold = 0.001; // must be adjusted according to the noise leve set by digitizer + fPRELocMaxCut = 0.03 ; + fECALocMaxCut = 0.03 ; + fHCALocMaxCut = 0.03 ; - fW0 = 4.5 ; - fW0CPV = 4.0 ; + fPREW0 = 4.0 ; + fECAW0 = 4.5 ; + fHCAW0 = 4.5 ; fTimeGate = 1.e-8 ; fToUnfold = kFALSE ; - - TString clusterizerName( GetName()) ; - if (clusterizerName.IsNull() ) - clusterizerName = "Default" ; - clusterizerName.Append(":") ; - clusterizerName.Append(Version()) ; - SetName(clusterizerName) ; - fRecPointsInRun = 0 ; - + + fRecPointsInRun = 0 ; + } //____________________________________________________________________________ @@ -371,7 +331,7 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster // which is compared to a digit (d2) not yet in a cluster - AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ; + AliEMCALGeometry * geom = AliEMCALGetter::Instance()->EMCALGeometry() ; Int_t rv = 0 ; @@ -380,8 +340,9 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d Int_t relid2[4] ; geom->AbsToRelNumbering(d2->GetId(), relid2) ; - - if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same EMCAL Arm + + if ( (relid1[0] == relid2[0]) && // inside the same EMCAL Arm + (relid1[1]==relid2[1]) ) { // and same tower section Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ; Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ; @@ -398,33 +359,25 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d else { if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) ) - rv=2 ; + rv=0 ; } - return rv ; -} -//____________________________________________________________________________ -Bool_t AliEMCALClusterizerv1::IsInTower(AliEMCALDigit * digit) const -{ - // Tells if (true) or not (false) the digit is in a EMCAL-Tower - - Bool_t rv = kFALSE ; - if (!digit->IsInPreShower()) - rv = kTRUE; + if (gDebug == 2 ) + Info("AreNeighbours", "neighbours=%d, id1=%d, relid1=%d,%d,%d,%d \n id2=%d, relid2=%d,%d,%d,%d", + rv, d1->GetId(), relid1[0], relid1[1], relid1[2], relid1[3], + d2->GetId(), relid2[0], relid2[1], relid2[2], relid2[3]) ; + return rv ; } //____________________________________________________________________________ -Bool_t AliEMCALClusterizerv1::IsInPreShower(AliEMCALDigit * digit) const +void AliEMCALClusterizerv1::Unload() { - // Tells if (true) or not (false) the digit is in a EMCAL-PreShower - - Bool_t rv = kFALSE ; - if (digit->IsInPreShower()) - rv = kTRUE; - return rv ; + AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + gime->EmcalLoader()->UnloadDigits() ; + gime->EmcalLoader()->UnloadRecPoints() ; } - + //____________________________________________________________________________ void AliEMCALClusterizerv1::WriteRecPoints(Int_t event) { @@ -432,74 +385,70 @@ void AliEMCALClusterizerv1::WriteRecPoints(Int_t event) // Creates new branches with given title // fills and writes into TreeR. - AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; - TObjArray * towerRecPoints = gime->TowerRecPoints() ; - TObjArray * preshoRecPoints = gime->PreShowerRecPoints() ; + AliEMCALGetter *gime = AliEMCALGetter::Instance() ; + + TObjArray * aPRERecPoints = gime->PRERecPoints() ; + TObjArray * aECARecPoints = gime->ECARecPoints() ; + TObjArray * aHCARecPoints = gime->HCARecPoints() ; + TClonesArray * digits = gime->Digits() ; - TTree * treeR ; + TTree * treeR = gime->TreeR(); ; - if(fToSplit){ - if(!fSplitFile) - return ; - fSplitFile->cd() ; - TString name("TreeR") ; - name += event ; - treeR = dynamic_cast(fSplitFile->Get(name)); - } - else{ - treeR = gAlice->TreeR(); - } - - if(!treeR){ - gAlice->MakeTree("R", fSplitFile); - treeR = gAlice->TreeR() ; - } - Int_t index ; - //Evaluate position, dispersion and other RecPoint properties... - for(index = 0; index < towerRecPoints->GetEntries(); index++) - (dynamic_cast(towerRecPoints->At(index)))->EvalAll(fW0,digits) ; + + //Evaluate position, dispersion and other RecPoint properties for PRE section + for(index = 0; index < aPRERecPoints->GetEntries(); index++) + (dynamic_cast(aPRERecPoints->At(index)))->EvalAll(fPREW0,digits) ; + aPRERecPoints->Sort() ; - towerRecPoints->Sort() ; + for(index = 0; index < aPRERecPoints->GetEntries(); index++) + (dynamic_cast(aPRERecPoints->At(index)))->SetIndexInList(index) ; + + aPRERecPoints->Expand(aPRERecPoints->GetEntriesFast()) ; + + //Evaluate position, dispersion and other RecPoint properties for EC section + for(index = 0; index < aECARecPoints->GetEntries(); index++) + (dynamic_cast(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ; + + aECARecPoints->Sort() ; - for(index = 0; index < towerRecPoints->GetEntries(); index++) - (dynamic_cast(towerRecPoints->At(index)))->SetIndexInList(index) ; + for(index = 0; index < aECARecPoints->GetEntries(); index++) + (dynamic_cast(aECARecPoints->At(index)))->SetIndexInList(index) ; - towerRecPoints->Expand(towerRecPoints->GetEntriesFast()) ; - //Now the same for pre shower - for(index = 0; index < preshoRecPoints->GetEntries(); index++) - (dynamic_cast(preshoRecPoints->At(index)))->EvalAll(fW0CPV,digits) ; - preshoRecPoints->Sort() ; + aECARecPoints->Expand(aECARecPoints->GetEntriesFast()) ; + + //Evaluate position, dispersion and other RecPoint properties for HCA section + for(index = 0; index < aHCARecPoints->GetEntries(); index++) + (dynamic_cast(aHCARecPoints->At(index)))->EvalAll(fHCAW0,digits) ; + + aHCARecPoints->Sort() ; - for(index = 0; index < preshoRecPoints->GetEntries(); index++) - (dynamic_cast(preshoRecPoints->At(index)))->SetIndexInList(index) ; + for(index = 0; index < aHCARecPoints->GetEntries(); index++) + (dynamic_cast(aHCARecPoints->At(index)))->SetIndexInList(index) ; - preshoRecPoints->Expand(preshoRecPoints->GetEntriesFast()) ; - + aHCARecPoints->Expand(aHCARecPoints->GetEntriesFast()) ; + Int_t bufferSize = 32000 ; Int_t splitlevel = 0 ; - - //First Tower branch - TBranch * towerBranch = treeR->Branch("EMCALTowerRP","TObjArray",&towerRecPoints,bufferSize,splitlevel); - towerBranch->SetTitle(BranchName()); - //Now Pre Shower branch - TBranch * preshoBranch = treeR->Branch("EMCALPreShoRP","TObjArray",&preshoRecPoints,bufferSize,splitlevel); - preshoBranch->SetTitle(BranchName()); - - //And Finally clusterizer branch - AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ; - TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1", - &cl,bufferSize,splitlevel); - clusterizerBranch->SetTitle(BranchName()); - - towerBranch ->Fill() ; - preshoBranch ->Fill() ; - clusterizerBranch->Fill() ; - - treeR->AutoSave() ; //Write(0,kOverwrite) ; - if(gAlice->TreeR()!=treeR) - treeR->Delete(); + //PRE section branch + TBranch * branchPRE = treeR->Branch("EMCALPRERP","TObjArray",&aPRERecPoints,bufferSize,splitlevel); + branchPRE->SetTitle(BranchName()); + + //EC section branch + TBranch * branchECA = treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel); + branchECA->SetTitle(BranchName()); + + //HCA section branch + TBranch * branchHCA = treeR->Branch("EMCALHCARP","TObjArray",&aHCARecPoints,bufferSize,splitlevel); + branchHCA->SetTitle(BranchName()); + + branchPRE->Fill() ; + branchECA->Fill() ; + branchHCA->Fill() ; + + gime->WriteRecPoints("OVERWRITE"); + gime->WriteClusterizer("OVERWRITE"); } //____________________________________________________________________________ @@ -508,106 +457,266 @@ void AliEMCALClusterizerv1::MakeClusters() // Steering method to construct the clusters stored in a list of Reconstructed Points // A cluster is defined as a list of neighbour digits - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; - - TObjArray * towerRecPoints = gime->TowerRecPoints(BranchName()) ; - TObjArray * preshoRecPoints = gime->PreShowerRecPoints(BranchName()) ; - towerRecPoints->Delete() ; - preshoRecPoints->Delete() ; - + AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + + AliEMCALGeometry * geom = gime->EMCALGeometry() ; + + + TObjArray * aPRERecPoints = gime->PRERecPoints() ; + TObjArray * aECARecPoints = gime->ECARecPoints() ; + TObjArray * aHCARecPoints = gime->HCARecPoints() ; + + aPRERecPoints->Delete() ; + aECARecPoints->Delete() ; + aHCARecPoints->Delete() ; + TClonesArray * digits = gime->Digits() ; - if ( !digits ) { - Fatal("MakeClusters -> Digits with name %s not found", GetName() ) ; - } + + TIter next(digits) ; + AliEMCALDigit * digit ; + Int_t ndigECA=0, ndigPRE=0, ndigHCA=0 ; + + // count the number of digits in ECA section + while ( (digit = dynamic_cast(next())) ) { // scan over the list of digits + if (geom->IsInECA(digit->GetId())) + ndigECA++ ; + else if (geom->IsInPRE(digit->GetId())) + ndigPRE++; + else if (geom->IsInHCA(digit->GetId())) + ndigHCA++; + else { + Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ; + abort() ; + } + } + + // add amplitude of PRE and ECA sections + Int_t digECA ; + for (digECA = 0 ; digECA < ndigECA ; digECA++) { + digit = dynamic_cast(digits->At(digECA)) ; + Int_t digPRE ; + for (digPRE = ndigECA ; digPRE < ndigECA+ndigPRE ; digPRE++) { + AliEMCALDigit * digitPRE = dynamic_cast(digits->At(digPRE)) ; + if ( geom->AreInSameTower(digit->GetId(), digitPRE->GetId()) ){ + Float_t amp = static_cast(digit->GetAmp()) + geom->GetSummationFraction() * static_cast(digitPRE->GetAmp()) + 0.5 ; + digit->SetAmp(static_cast(amp)) ; + break ; + } + } + if (gDebug) + Info("MakeClusters","id = %d amp = %d", digit->GetId(), digit->GetAmp()) ; + } + TClonesArray * digitsC = dynamic_cast(digits->Clone()) ; // Clusterization starts TIter nextdigit(digitsC) ; - AliEMCALDigit * digit ; - Bool_t notremoved = kTRUE ; + Bool_t notremovedECA = kTRUE, notremovedPRE = kTRUE ; while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC AliEMCALRecPoint * clu = 0 ; - TArrayI clusterdigitslist(1500) ; - Int_t index ; - if (( IsInTower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fTowerClusteringThreshold ) || - ( IsInPreShower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fPreShoClusteringThreshold ) ) { - - Int_t iDigitInCluster = 0 ; + TArrayI clusterPREdigitslist(50), clusterECAdigitslist(50), clusterHCAdigitslist(50); + + Bool_t inPRE = kFALSE, inECA = kFALSE, inHCA = kFALSE ; + if( geom->IsInPRE(digit->GetId()) ) { + inPRE = kTRUE ; + } + else if( geom->IsInECA(digit->GetId()) ) { + inECA = kTRUE ; + } + else if( geom->IsInHCA(digit->GetId()) ) { + inHCA = kTRUE ; + } + + if (gDebug == 2) { + if (inPRE) + Info("MakeClusters","id = %d, ene = %f , thre = %f ", + digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ; + if (inECA) + Info("MakeClusters","id = %d, ene = %f , thre = %f", + digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ; + if (inHCA) + Info("MakeClusters","id = %d, ene = %f , thre = %f", + digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold ) ; + } + + if ( (inPRE && (Calibrate(digit->GetAmp(), 0) > fPREClusteringThreshold )) || + (inECA && (Calibrate(digit->GetAmp(), 1) > fECAClusteringThreshold )) || + (inHCA && (Calibrate(digit->GetAmp(), 2) > fHCAClusteringThreshold )) ) { - if ( IsInTower(digit) ) { + Int_t iDigitInPRECluster = 0, iDigitInECACluster = 0, iDigitInHCACluster = 0; + Int_t where ; // PRE = 0, ECAl = 1, HCAL = 2 + + // Find the seed in each of the section ECAL/PRE/HCAL + + if( geom->IsInECA(digit->GetId()) ) { + where = 1 ; // to tell we are in ECAL // start a new Tower RecPoint - if(fNumberOfTowerClusters >= towerRecPoints->GetSize()) - towerRecPoints->Expand(2*fNumberOfTowerClusters+1) ; - - towerRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfTowerClusters) ; - clu = dynamic_cast(towerRecPoints->At(fNumberOfTowerClusters)) ; - fNumberOfTowerClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower())) ; - clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; - iDigitInCluster++ ; + if(fNumberOfECAClusters >= aECARecPoints->GetSize()) + aECARecPoints->Expand(2*fNumberOfECAClusters+1) ; + AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ; + rp->SetECA() ; + aECARecPoints->AddAt(rp, fNumberOfECAClusters) ; + clu = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; + fNumberOfECAClusters++ ; + clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)) ; + clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ; + iDigitInECACluster++ ; digitsC->Remove(digit) ; + if (gDebug == 2 ) + Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ; - } else { - + } + else if( geom->IsInPRE(digit->GetId()) ) { + where = 0 ; // to tell we are in PRE // start a new Pre Shower cluster - if(fNumberOfPreShoClusters >= preshoRecPoints->GetSize()) - preshoRecPoints->Expand(2*fNumberOfPreShoClusters+1); - - preshoRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfPreShoClusters) ; + if(fNumberOfPREClusters >= aPRERecPoints->GetSize()) + aPRERecPoints->Expand(2*fNumberOfPREClusters+1); + AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ; + rp->SetPRE() ; + aPRERecPoints->AddAt(rp, fNumberOfPREClusters) ; + clu = dynamic_cast(aPRERecPoints->At(fNumberOfPREClusters)) ; + fNumberOfPREClusters++ ; + clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)); + clusterPREdigitslist[iDigitInPRECluster] = digit->GetIndexInList() ; + iDigitInPRECluster++ ; + digitsC->Remove(digit) ; + if (gDebug == 2 ) + Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ; - clu = dynamic_cast(preshoRecPoints->At(fNumberOfPreShoClusters)) ; - fNumberOfPreShoClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower() ) ); - clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; - iDigitInCluster++ ; - digitsC->Remove(digit) ; nextdigit.Reset() ; + + // Here we remove remaining ECA digits, which cannot make a cluster - // Here we remove remaining Tower digits, which cannot make a cluster - - if( notremoved ) { + if( notremovedECA ) { while( ( digit = dynamic_cast(nextdigit()) ) ) { - if( IsInTower(digit) ) + if( geom->IsInECA(digit->GetId()) ) digitsC->Remove(digit) ; else break ; } - notremoved = kFALSE ; + notremovedECA = kFALSE ; } + + } + else if( geom->IsInHCA(digit->GetId()) ) { + where = 2 ; // to tell we are in HCAL + // start a new HCAL cluster + if(fNumberOfHCAClusters >= aHCARecPoints->GetSize()) + aHCARecPoints->Expand(2*fNumberOfHCAClusters+1); + AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ; + rp->SetHCA() ; + aHCARecPoints->AddAt(rp, fNumberOfHCAClusters) ; + clu = dynamic_cast(aHCARecPoints->At(fNumberOfHCAClusters)) ; + fNumberOfHCAClusters++ ; + clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)); + clusterHCAdigitslist[iDigitInHCACluster] = digit->GetIndexInList() ; + iDigitInHCACluster++ ; + digitsC->Remove(digit) ; + if (gDebug == 2 ) + Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold) ; + + nextdigit.Reset() ; + + // Here we remove remaining PRE digits, which cannot make a cluster - } // else + if( notremovedPRE ) { + while( ( digit = dynamic_cast(nextdigit()) ) ) { + if( geom->IsInPRE(digit->GetId()) ) + digitsC->Remove(digit) ; + else + break ; + } + notremovedPRE = kFALSE ; + } + } nextdigit.Reset() ; AliEMCALDigit * digitN ; - index = 0 ; - while (index < iDigitInCluster){ // scan over digits already in cluster - digit = (AliEMCALDigit*)digits->At(clusterdigitslist[index]) ; + Int_t index = 0 ; + + // Do the Clustering in each of the three section ECAL/PRE/HCAL + + while (index < iDigitInECACluster){ // scan over digits already in cluster + digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]) ; index++ ; while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! + // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ; switch (ineb ) { case 0 : // not a neighbour break ; case 1 : // are neighbours - clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->IsInPreShower() ) ) ; - clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; - iDigitInCluster++ ; + clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 1) ) ; + clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; + iDigitInECACluster++ ; digitsC->Remove(digitN) ; break ; case 2 : // too far from each other - goto endofloop; + goto endofloop1; } // switch } // while digitN - endofloop: ; + endofloop1: ; nextdigit.Reset() ; - } // loop over cluster + } // loop over ECA cluster + + index = 0 ; + while (index < iDigitInPRECluster){ // scan over digits already in cluster + digit = (AliEMCALDigit*)digits->At(clusterPREdigitslist[index]) ; + index++ ; + while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits + Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! + // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ; + switch (ineb ) { + case 0 : // not a neighbour + break ; + case 1 : // are neighbours + clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 0) ) ; + clusterPREdigitslist[iDigitInPRECluster] = digitN->GetIndexInList() ; + iDigitInPRECluster++ ; + digitsC->Remove(digitN) ; + break ; + case 2 : // too far from each other + goto endofloop2; + } // switch + + } // while digitN + + endofloop2: ; + nextdigit.Reset() ; + } // loop over PRE cluster + + index = 0 ; + while (index < iDigitInHCACluster){ // scan over digits already in cluster + digit = (AliEMCALDigit*)digits->At(clusterHCAdigitslist[index]) ; + index++ ; + while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits + Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! + //Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ; + switch (ineb ) { + case 0 : // not a neighbour + break ; + case 1 : // are neighbours + clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 2) ) ; + clusterHCAdigitslist[iDigitInHCACluster] = digitN->GetIndexInList() ; + iDigitInHCACluster++ ; + digitsC->Remove(digitN) ; + break ; + case 2 : // too far from each other + goto endofloop3; + } // switch + } // while digitN + + endofloop3: ; + nextdigit.Reset() ; + } // loop over HCA cluster + } // energy theshold } // while digit delete digitsC ; @@ -675,18 +784,24 @@ void AliEMCALClusterizerv1::Print(Option_t * option)const message += taskName.Data() ; message += "\n Branch: " ; message += GetName() ; - message += "\n EMC Clustering threshold = " ; - message += fTowerClusteringThreshold ; - message += "\n EMC Local Maximum cut = " ; - message += fTowerLocMaxCut ; - message += "\n EMC Logarothmic weight = " ; - message += fW0 ; - message += "\n CPV Clustering threshold = " ; - message += fPreShoClusteringThreshold ; - message += "\n CPV Local Maximum cut = " ; - message += fPreShoLocMaxCut ; - message += "\n CPV Logarothmic weight = " ; - message += fW0CPV ; + message += "\n Pre Shower Clustering threshold = " ; + message += fPREClusteringThreshold ; + message += "\n Pre Shower Local Maximum cut = " ; + message += fPRELocMaxCut ; + message += "\n Pre Shower Logarothmic weight = " ; + message += fPREW0 ; + message += "\n ECA Clustering threshold = " ; + message += fECAClusteringThreshold ; + message += "\n ECA Local Maximum cut = " ; + message += fECALocMaxCut ; + message += "\n ECA Logarothmic weight = " ; + message += fECAW0 ; + message += "\n Pre Shower Clustering threshold = " ; + message += fHCAClusteringThreshold ; + message += "\n HCA Local Maximum cut = " ; + message += fHCALocMaxCut ; + message += "\n HCA Logarothmic weight = " ; + message += fHCAW0 ; if(fToUnfold) message +="\nUnfolding on\n" ; else @@ -705,80 +820,98 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) { // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer - TObjArray * towerRecPoints = AliEMCALGetter::GetInstance()->TowerRecPoints() ; - TObjArray * preshoRecPoints = AliEMCALGetter::GetInstance()->PreShowerRecPoints() ; - - TString message("\n") ; + TObjArray * aPRERecPoints = AliEMCALGetter::Instance()->PRERecPoints() ; + TObjArray * aECARecPoints = AliEMCALGetter::Instance()->ECARecPoints() ; + TObjArray * aHCARecPoints = AliEMCALGetter::Instance()->HCARecPoints() ; - message += "event " ; - message += gAlice->GetEvNumber() ; - message += "\n Found " ; - message += towerRecPoints->GetEntriesFast() ; - message += " TOWER Rec Points and " ; - message += preshoRecPoints->GetEntriesFast() ; - message += " PRE SHOWER RecPoints\n" ; - - fRecPointsInRun += towerRecPoints->GetEntriesFast() ; - fRecPointsInRun += preshoRecPoints->GetEntriesFast() ; + Info("PrintRecPoints", "Clusterization result:") ; + + printf("event # %d\n", gAlice->GetEvNumber() ) ; + printf(" Found %d PRE SHOWER RecPoints, %d ECA Rec Points and %d HCA Rec Points\n ", + aPRERecPoints->GetEntriesFast(), aECARecPoints->GetEntriesFast(), aHCARecPoints->GetEntriesFast() ) ; - char * tempo = new char[8192]; + fRecPointsInRun += aPRERecPoints->GetEntriesFast() ; + fRecPointsInRun += aECARecPoints->GetEntriesFast() ; + fRecPointsInRun += aHCARecPoints->GetEntriesFast() ; if(strstr(option,"all")) { - message += "Tower clusters\n" ; - message += "Index Ene(MeV) Multi Module phi r theta Lambda 1 Lambda 2 # of prim Primaries list\n" ; - + //Pre shower recPoints + + printf("-----------------------------------------------------------------------\n") ; + printf("Clusters in PRE section\n") ; + printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + Int_t index ; - for (index = 0 ; index < towerRecPoints->GetEntries() ; index++) { - AliEMCALTowerRecPoint * rp = dynamic_cast(towerRecPoints->At(index)) ; + + for (index = 0 ; index < aPRERecPoints->GetEntries() ; index++) { + AliEMCALTowerRecPoint * rp = dynamic_cast(aPRERecPoints->At(index)) ; TVector3 globalpos; rp->GetGlobalPosition(globalpos); + TVector3 localpos; + rp->GetLocalPosition(localpos); + Float_t lambda[2]; + rp->GetElipsAxis(lambda); + Int_t * primaries; + Int_t nprimaries; + primaries = rp->GetPrimaries(nprimaries); + printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", + rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), + globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), + rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; + for (Int_t iprimary=0; iprimaryGetEntries() ; index++) { + AliEMCALTowerRecPoint * rp = dynamic_cast(aECARecPoints->At(index)) ; + TVector3 globalpos; + rp->GetGlobalPosition(globalpos); + TVector3 localpos; + rp->GetLocalPosition(localpos); Float_t lambda[2]; rp->GetElipsAxis(lambda); Int_t * primaries; Int_t nprimaries; primaries = rp->GetPrimaries(nprimaries); - sprintf(tempo, "\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ", - rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), - globalpos.X(), globalpos.Y(), globalpos.Z(), lambda[0], lambda[1], nprimaries) ; - message += tempo ; - + printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", + rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), + globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), + rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; for (Int_t iprimary=0; iprimaryGetEntries() ; index++) { - AliEMCALTowerRecPoint * rp = dynamic_cast(preshoRecPoints->At(index)) ; + for (index = 0 ; index < aHCARecPoints->GetEntries() ; index++) { + AliEMCALTowerRecPoint * rp = dynamic_cast(aHCARecPoints->At(index)) ; TVector3 globalpos; rp->GetGlobalPosition(globalpos); + TVector3 localpos; + rp->GetLocalPosition(localpos); Float_t lambda[2]; rp->GetElipsAxis(lambda); - Int_t * primaries; + Int_t * primaries; Int_t nprimaries; primaries = rp->GetPrimaries(nprimaries); - sprintf(tempo, "\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ", - rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), - globalpos.X(), globalpos.Y(), globalpos.Z(), lambda[0], lambda[1], nprimaries) ; - message += tempo ; - + printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", + rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), + globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), + rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; for (Int_t iprimary=0; iprimary