[u/mrichter/AliRoot.git] / MUON / TestRecPoints.C

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

/* $Id: TestRecPoints.C 23207 2007-12-20 09:59:20Z ivana $ */

/// \ingroup macros
/// \file TestRecPoints.C
/// \brief This macro is to be used to check the trigger and tracking chamber behaviour
/// through the RecPoints analysis
///
/// \author D. Stocco (Torino)

#if !defined(__CINT__) || defined(__MAKECINT__)
// ROOT
#include "TString.h"
#include "TH1.h"
#include "TList.h"
#include "TFile.h"
#include "TDirectory.h"

#include "TTree.h"
#include "TClonesArray.h"

// STEER
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliLoader.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliCDBPath.h"

// tracker
#include "AliGeomManager.h"

// MUON
#include "AliMpCDB.h"
#include "AliMpDDL.h"
#include "AliMpDDLStore.h"
#include "AliMpConstants.h"

// trigger
#include "AliMUONVTriggerStore.h"
#include "AliMUONDigitStoreV1.h"
#include "AliMUONDigitMaker.h"
#include "AliMUONLocalTrigger.h"

// tracker
#include "AliMUONClusterStoreV2.h"
#include "AliMUONClusterFinderMLEM.h"
#include "AliMUONPreClusterFinder.h"
#include "AliMUONSimpleClusterServer.h"
#include "AliMUONGeometryTransformer.h"
#include "AliMUONVDigitStore.h"
#include "AliMUONDigitStoreV2R.h"
#include "AliMUONVCluster.h"
#include "AliMUONRecoParam.h"

#endif

const Int_t kNtrigChambers = AliMpConstants::NofTriggerChambers();
const Int_t kNcathodes = AliMpConstants::NofCathodes();

const Int_t kNtrackChambers = AliMpConstants::NofTrackingChambers();

enum {kOnlyTracker, kOnlyTrigger, kTrackTrig};

Int_t GetPlane(Int_t ch, Int_t cath){return kNcathodes * ch + cath;}
void ClusterSize(TList&, AliMUONVDigit*, Int_t&, Int_t);

// Main Method
void TestRecPoints(TString baseDir=".", TString outDir=".", Float_t adcCut = 10., Int_t whatToTest=kTrackTrig, Int_t runNumber=0, TString cdbStorage="local://$ALICE_ROOT")
{
  const Int_t kNplanes = kNtrigChambers * kNcathodes;
  const Int_t kNslats = 18;
  
  TList histoListTrig;
  
  TH1F* trigStripMult[kNplanes];
  TH1F* trigClusterSize[kNplanes];
  TH1F* trigClusterMult[kNplanes];
  TString histoName, histoTitle;
  
  if(whatToTest!=kOnlyTracker){
    for(Int_t ch=0; ch<kNtrigChambers; ch++){
      for(Int_t cath=0; cath<kNcathodes; cath++){
        Int_t iplane = GetPlane(ch, cath);
        histoName = "stripMult";
        histoName.Append(Form("Ch%iCath%i",ch,cath));
        histoTitle = Form("Strip multiplicity Ch %i Cathode %i", ch, cath);
        trigStripMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
        trigStripMult[iplane]->SetXTitle("slat");
        histoListTrig.Add(trigStripMult[iplane]);
        
        histoName = "clusterSize";
        histoName.Append(Form("Ch%iCath%i",ch,cath));
        histoTitle = Form("Cluster size Ch %i Cathode %i", ch, cath);
        trigClusterSize[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), 10, 0.-0.5, (Float_t)10 - 0.5);
        trigClusterSize[iplane]->SetXTitle("cluster size");
        histoListTrig.Add(trigClusterSize[iplane]);
        
        histoName = "clusterMult";
        histoName.Append(Form("Ch%iCath%i",ch,cath));
        histoTitle = Form("Cluster multiplicity Ch %i Cathode %i", ch, cath);
        trigClusterMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
        trigClusterMult[iplane]->SetXTitle("slat");
        histoListTrig.Add(trigClusterMult[iplane]);
      }
    }
  }
  
  const Int_t kMaxDetElem = 26;
  TList histoListTrack;
  TH1F* trackClusterSize[kNtrackChambers];
  TH1F* trackClusterMult[kNtrackChambers];
  TH1F* trackClusterYield[kNtrackChambers];
  
  if(whatToTest!=kOnlyTrigger){
    for(Int_t ch=0; ch<kNtrackChambers; ch++){
      histoName = "clusterSize";
      histoName.Append(Form("Ch%i",ch));
      histoTitle = Form("Cluster size Ch %i", ch);
      trackClusterSize[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
      trackClusterSize[ch]->SetXTitle("cluster size");
      histoListTrack.Add(trackClusterSize[ch]);
      
      histoName = "clusterMult";
      histoName.Append(Form("Ch%i",ch));
      histoTitle = Form("Cluster multiplicity vs. Charge Ch %i", ch);
      trackClusterMult[ch] = new TH1F(histoName.Data(), histoTitle.Data(), 200, 0, 15000.);
      trackClusterMult[ch]->SetXTitle("Charge (ADC)");
      histoListTrack.Add(trackClusterMult[ch]);
      
      histoName = "clusterYield";
      histoName.Append(Form("Ch%i",ch));
      histoTitle = Form("Cluster yield vs. DetElem Ch %i", ch);
      trackClusterYield[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
      trackClusterYield[ch]->SetXTitle("Detector element #");
      histoListTrack.Add(trackClusterYield[ch]);
    }
  }
  
  
  // Creating Run Loader and opening RecPoints
  TString filename = baseDir.Data();
  filename.Append("/galice.root");
  AliRunLoader * RunLoader = AliRunLoader::Open(filename.Data(),"MUONLoader","UPDATE");
  if (RunLoader ==0x0) {
    printf(">>> Error : Error Opening %s file \n",filename.Data());
    return;
  }
  
  // Loading MUON subsystem
  AliLoader* MUONLoader = RunLoader->GetDetectorLoader("MUON");
  if(whatToTest!=kOnlyTracker)  MUONLoader->LoadRecPoints("READ");
  if(whatToTest!=kOnlyTrigger)  MUONLoader->LoadDigits("READ");
  
  TTree *treeR = 0x0, *treeD = 0x0;
  AliMUONVTriggerStore* trigStore = 0x0;
  AliMUONLocalTrigger* locTrg = 0x0;
  
  AliMUONVDigitStore* digitStoreTrack = 0x0;
  AliMUONDigitStoreV2R digitStoreTrackCut;
  AliMUONVCluster* cluster = 0x0;
  
  // Load segmentation
  AliCDBManager::Instance()->SetDefaultStorage(cdbStorage.Data());
  AliCDBManager::Instance()->SetRun(runNumber);
  
  AliMpCDB::LoadDDLStore();
  
  AliMUONGeometryTransformer* transformer = 0x0;
  
  AliMUONRecoParam* recoParam = 0x0;
  
  AliMUONClusterStoreV2 clusterStore;
  AliMUONVClusterFinder* clusterFinder = 0x0;
  
  AliMUONSimpleClusterServer* clusterServer = 0x0;
  
  Int_t firstChamber(0);
  Int_t lastChamber(9);
  
  if(whatToTest!=kOnlyTrigger){
    // Import TGeo geometry
    TString geoFilename = baseDir.Data();
    geoFilename.Append("/geometry.root");
    if ( ! AliGeomManager::GetGeometry() ) {
      AliGeomManager::LoadGeometry(geoFilename);
      if (! AliGeomManager::GetGeometry() ) {
        printf(">>> Error : getting geometry from file %s failed\n", geoFilename.Data());
        return;
      }
    }
    transformer = new AliMUONGeometryTransformer();
    // Load geometry data
    transformer->LoadGeometryData();
    // Load reconstruction parameters
    AliCDBPath path("MUON","Calib","RecoParam");
    AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
    if(entry) {
      recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
      entry->SetOwner(0);
      AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
    }
    if (!recoParam) {
      printf("Couldn't find RecoParam object in OCDB: create default one");
      recoParam = AliMUONRecoParam::GetLowFluxParam();
    }
    recoParam->Print("FULL");
    clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
    clusterServer = new AliMUONSimpleClusterServer(clusterFinder,*transformer);
  }
  
  AliMUONDigitStoreV1 digitStore;
  AliMUONVDigit* mDigit;
  
  Int_t clusterSize;
  
  AliMUONDigitMaker digitMaker;
  TList digitsList[kNplanes];
  
  Int_t nevents = RunLoader->GetNumberOfEvents();
  Int_t analysisFrac = nevents/10+1;
  
  printf("\nNumber of events = %i\n\n",nevents);
  
  for(Int_t ievent=0; ievent<nevents; ievent++){
    if(ievent%analysisFrac==0) printf("Analysing event = %i\n",ievent);
    RunLoader->GetEvent(ievent);
    if(whatToTest!=kOnlyTracker){
      digitStore.Clear();
      treeR = MUONLoader->TreeR();
      trigStore = AliMUONVTriggerStore::Create(*treeR);
      
      if ( trigStore == 0x0 ) continue;
      trigStore->Clear();
      trigStore->Connect(*treeR);
      treeR->GetEvent(0);
      
      TIter nextLocal(trigStore->CreateLocalIterator());
      while ( (locTrg = static_cast<AliMUONLocalTrigger*>( nextLocal() )) != NULL )
      {
        TArrayS xyPattern[2];
        locTrg->GetXPattern(xyPattern[0]);
        locTrg->GetYPattern(xyPattern[1]);
        
        Int_t nBoard = locTrg->LoCircuit();
        digitMaker.TriggerDigits(nBoard, xyPattern, digitStore);
      }
      
      TIter next(digitStore.CreateIterator());
      
      while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
      {
        Int_t detElemId = mDigit->DetElemId();
        Int_t ch = detElemId/100 - 11;
        Int_t cathode = mDigit->Cathode();
        Int_t slat = detElemId%100;
        Int_t iplane = GetPlane(ch, cathode);
        trigStripMult[iplane]->Fill(slat);
        digitsList[iplane].Add(mDigit);
      } // loop on digits
      
      for(Int_t iplane=0; iplane<kNplanes; iplane++){
        while(digitsList[iplane].GetEntries()){
          clusterSize=1;
          mDigit = (AliMUONVDigit*)digitsList[iplane].At(0);
          digitsList[iplane].Remove(mDigit);
          ClusterSize(digitsList[iplane],mDigit,clusterSize,1);
          //if(clusterSize>1) printf("Cluster size = %i\n\n",clusterSize);
          trigClusterSize[iplane]->Fill(clusterSize);
          
          Int_t detElemId = mDigit->DetElemId();
          Int_t slat = detElemId%100;
          trigClusterMult[iplane]->Fill(slat);
        } // loop o sorted digits
      } // loop on planes
    } // trigger part
    
    if(whatToTest!=kOnlyTrigger){
      clusterStore.Clear();
      treeD = MUONLoader->TreeD();
      digitStoreTrack = AliMUONVDigitStore::Create(*treeD);
      
      if ( digitStoreTrack == 0x0 ) continue;
      digitStoreTrack->Clear();
      digitStoreTrackCut.Clear();
      digitStoreTrack->Connect(*treeD);
      treeD->GetEvent(0);
      
      // Cut low charge channels (pedestal subtraction)
      TIter nextDigitTrack(digitStoreTrack->CreateIterator());
      while ( ( mDigit = static_cast<AliMUONVDigit*>(nextDigitTrack()) ) )
      {
        //printf("Digit class = %s",mDigit->ClassName());
        Float_t charge = mDigit->Charge();
        if(charge<adcCut) continue;
        digitStoreTrackCut.Add(*mDigit, AliMUONVDigitStore::kDeny);
      } // loop on digits
      
      TIter nextDigitTrackCut(digitStoreTrackCut.CreateIterator());
      clusterServer->UseDigits(nextDigitTrackCut);
      
      for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
        clusterServer->Clusterize(ch, clusterStore, AliMpArea(),recoParam);
      }
      
      TIter next(clusterStore.CreateIterator());
      while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) )
      {
        Float_t charge = cluster->GetCharge();
        if(charge==0.) continue;
        Int_t ch = cluster->GetChamberId();
        Int_t npads = cluster->GetNDigits();
        Int_t detElemId = cluster->GetDetElemId()%100;
        //printf("charge = %f   pads = %i   detElemId = %i\n",charge,npads,detElemId);
        
        trackClusterSize[ch]->Fill(npads);
        trackClusterMult[ch]->Fill(charge);
        trackClusterYield[ch]->Fill(detElemId);      
      } // loop on clusters
    } // tracker part
  } // loop on events
  
  MUONLoader->UnloadRecPoints();
  MUONLoader->UnloadDigits();
  RunLoader->UnloadAll();
  delete RunLoader;
  TString outFileName = outDir.Data();
  outFileName.Append("/outTestRecPoints.root");
  TFile* outFile = new TFile(outFileName.Data(), "RECREATE");
  TDirectory* dir = 0x0;
  if(whatToTest!=kOnlyTracker){
    outFile->cd();
    dir = outFile->mkdir("Trigger");
    dir->cd();
    histoListTrig.Write();
  }
  if(whatToTest!=kOnlyTrigger){
    outFile->cd();
    dir = outFile->mkdir("Tracker");
    dir->cd();
    histoListTrack.Write();
  }
  outFile->Close();
  printf("\nSee results in %s\n",outFileName.Data());
}

void ClusterSize(TList& list, AliMUONVDigit* refDigit, Int_t& clusterSize, Int_t depthLevel)
{
  AliMUONVDigit* mDigit = 0x0;
  for(Int_t idigit=0; idigit<list.GetEntries(); idigit++){
    mDigit = (AliMUONVDigit*) list.At(idigit);
    if(mDigit->DetElemId() != refDigit->DetElemId()) continue;
    Int_t diffX = TMath::Abs(mDigit->PadX() - refDigit->PadX());
    Int_t diffY = TMath::Abs(mDigit->PadY() - refDigit->PadY());
    if(diffX>1) continue;
    if(diffY>1) continue;
    if(diffX*diffY != 0) continue;
    clusterSize++;
    list.Remove(mDigit);
    //printf("DetElemId = %i   Level = %i   Ref. (%2i,%2i)    Found (%2i,%2i)\n",mDigit->DetElemId(),depthLevel,refDigit->PadX(),refDigit->PadY(),mDigit->PadX(),mDigit->PadY());
    ClusterSize(list, mDigit, clusterSize, depthLevel+1);
    Int_t val = idigit + depthLevel - clusterSize;
    //printf("Level = %i   Current digit = %i\t",depthLevel,idigit);
    idigit = TMath::Max(-1,val);
    //printf("New digit = %i\n",idigit);
  }
}
Commit	Line	Data
cc727848	1	/**************************************************************************
6f6a3f3b	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
cc727848	15
	16	/* $Id: TestRecPoints.C 23207 2007-12-20 09:59:20Z ivana $ */
	17
	18	/// \ingroup macros
	19	/// \file TestRecPoints.C
	20	/// \brief This macro is to be used to check the trigger and tracking chamber behaviour
	21	/// through the RecPoints analysis
	22	///
	23	/// \author D. Stocco (Torino)
	24
	25	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	26	// ROOT
	27	#include "TString.h"
	28	#include "TH1.h"
	29	#include "TList.h"
	30	#include "TFile.h"
	31	#include "TDirectory.h"
	32
	33	#include "TTree.h"
	34	#include "TClonesArray.h"
	35
	36	// STEER
	37	#include "AliRun.h"
	38	#include "AliRunLoader.h"
	39	#include "AliLoader.h"
	40	#include "AliCDBManager.h"
35be7ed7	41	#include "AliCDBEntry.h"
35be7ed7	42	#include "AliCDBPath.h"
cc727848	43
	44	// tracker
	45	#include "AliGeomManager.h"
	46
	47	// MUON
	48	#include "AliMpCDB.h"
	49	#include "AliMpDDL.h"
	50	#include "AliMpDDLStore.h"
	51	#include "AliMpConstants.h"
	52
	53	// trigger
	54	#include "AliMUONVTriggerStore.h"
	55	#include "AliMUONDigitStoreV1.h"
	56	#include "AliMUONDigitMaker.h"
	57	#include "AliMUONLocalTrigger.h"
	58
	59	// tracker
	60	#include "AliMUONClusterStoreV2.h"
	61	#include "AliMUONClusterFinderMLEM.h"
	62	#include "AliMUONPreClusterFinder.h"
	63	#include "AliMUONSimpleClusterServer.h"
	64	#include "AliMUONGeometryTransformer.h"
	65	#include "AliMUONVDigitStore.h"
6f6a3f3b	66	#include "AliMUONDigitStoreV2R.h"
cc727848	67	#include "AliMUONVCluster.h"
35be7ed7	68	#include "AliMUONRecoParam.h"
cc727848	69
	70	#endif
	71
	72	const Int_t kNtrigChambers = AliMpConstants::NofTriggerChambers();
	73	const Int_t kNcathodes = AliMpConstants::NofCathodes();
	74
	75	const Int_t kNtrackChambers = AliMpConstants::NofTrackingChambers();
	76
	77	enum {kOnlyTracker, kOnlyTrigger, kTrackTrig};
	78
	79	Int_t GetPlane(Int_t ch, Int_t cath){return kNcathodes * ch + cath;}
	80	void ClusterSize(TList&, AliMUONVDigit*, Int_t&, Int_t);
	81
	82	// Main Method
6f6a3f3b	83	void TestRecPoints(TString baseDir=".", TString outDir=".", Float_t adcCut = 10., Int_t whatToTest=kTrackTrig, Int_t runNumber=0, TString cdbStorage="local://$ALICE_ROOT")
cc727848	84	{
	85	const Int_t kNplanes = kNtrigChambers * kNcathodes;
	86	const Int_t kNslats = 18;
6f6a3f3b	87
cc727848	88	TList histoListTrig;
6f6a3f3b	89
cc727848	90	TH1F* trigStripMult[kNplanes];
	91	TH1F* trigClusterSize[kNplanes];
	92	TH1F* trigClusterMult[kNplanes];
	93	TString histoName, histoTitle;
6f6a3f3b	94
cc727848	95	if(whatToTest!=kOnlyTracker){
	96	for(Int_t ch=0; ch<kNtrigChambers; ch++){
	97	for(Int_t cath=0; cath<kNcathodes; cath++){
6f6a3f3b	98	Int_t iplane = GetPlane(ch, cath);
	99	histoName = "stripMult";
	100	histoName.Append(Form("Ch%iCath%i",ch,cath));
	101	histoTitle = Form("Strip multiplicity Ch %i Cathode %i", ch, cath);
	102	trigStripMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
	103	trigStripMult[iplane]->SetXTitle("slat");
	104	histoListTrig.Add(trigStripMult[iplane]);
	105
	106	histoName = "clusterSize";
	107	histoName.Append(Form("Ch%iCath%i",ch,cath));
	108	histoTitle = Form("Cluster size Ch %i Cathode %i", ch, cath);
	109	trigClusterSize[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), 10, 0.-0.5, (Float_t)10 - 0.5);
	110	trigClusterSize[iplane]->SetXTitle("cluster size");
	111	histoListTrig.Add(trigClusterSize[iplane]);
	112
	113	histoName = "clusterMult";
	114	histoName.Append(Form("Ch%iCath%i",ch,cath));
	115	histoTitle = Form("Cluster multiplicity Ch %i Cathode %i", ch, cath);
	116	trigClusterMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
	117	trigClusterMult[iplane]->SetXTitle("slat");
	118	histoListTrig.Add(trigClusterMult[iplane]);
cc727848	119	}
	120	}
	121	}
6f6a3f3b	122
cc727848	123	const Int_t kMaxDetElem = 26;
	124	TList histoListTrack;
	125	TH1F* trackClusterSize[kNtrackChambers];
	126	TH1F* trackClusterMult[kNtrackChambers];
	127	TH1F* trackClusterYield[kNtrackChambers];
6f6a3f3b	128
cc727848	129	if(whatToTest!=kOnlyTrigger){
	130	for(Int_t ch=0; ch<kNtrackChambers; ch++){
	131	histoName = "clusterSize";
	132	histoName.Append(Form("Ch%i",ch));
	133	histoTitle = Form("Cluster size Ch %i", ch);
	134	trackClusterSize[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
	135	trackClusterSize[ch]->SetXTitle("cluster size");
	136	histoListTrack.Add(trackClusterSize[ch]);
6f6a3f3b	137
cc727848	138	histoName = "clusterMult";
	139	histoName.Append(Form("Ch%i",ch));
	140	histoTitle = Form("Cluster multiplicity vs. Charge Ch %i", ch);
	141	trackClusterMult[ch] = new TH1F(histoName.Data(), histoTitle.Data(), 200, 0, 15000.);
6f6a3f3b	142	trackClusterMult[ch]->SetXTitle("Charge (ADC)");
cc727848	143	histoListTrack.Add(trackClusterMult[ch]);
6f6a3f3b	144
cc727848	145	histoName = "clusterYield";
	146	histoName.Append(Form("Ch%i",ch));
	147	histoTitle = Form("Cluster yield vs. DetElem Ch %i", ch);
	148	trackClusterYield[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
	149	trackClusterYield[ch]->SetXTitle("Detector element #");
	150	histoListTrack.Add(trackClusterYield[ch]);
	151	}
	152	}
	153
6f6a3f3b	154
cc727848	155	// Creating Run Loader and opening RecPoints
	156	TString filename = baseDir.Data();
	157	filename.Append("/galice.root");
	158	AliRunLoader * RunLoader = AliRunLoader::Open(filename.Data(),"MUONLoader","UPDATE");
	159	if (RunLoader ==0x0) {
	160	printf(">>> Error : Error Opening %s file \n",filename.Data());
	161	return;
	162	}
6f6a3f3b	163
cc727848	164	// Loading MUON subsystem
	165	AliLoader* MUONLoader = RunLoader->GetDetectorLoader("MUON");
	166	if(whatToTest!=kOnlyTracker) MUONLoader->LoadRecPoints("READ");
	167	if(whatToTest!=kOnlyTrigger) MUONLoader->LoadDigits("READ");
6f6a3f3b	168
cc727848	169	TTree treeR = 0x0, treeD = 0x0;
	170	AliMUONVTriggerStore* trigStore = 0x0;
	171	AliMUONLocalTrigger* locTrg = 0x0;
6f6a3f3b	172
cc727848	173	AliMUONVDigitStore* digitStoreTrack = 0x0;
6f6a3f3b	174	AliMUONDigitStoreV2R digitStoreTrackCut;
cc727848	175	AliMUONVCluster* cluster = 0x0;
6f6a3f3b	176
cc727848	177	// Load segmentation
	178	AliCDBManager::Instance()->SetDefaultStorage(cdbStorage.Data());
	179	AliCDBManager::Instance()->SetRun(runNumber);
6f6a3f3b	180
cc727848	181	AliMpCDB::LoadDDLStore();
6f6a3f3b	182
cc727848	183	AliMUONGeometryTransformer* transformer = 0x0;
6f6a3f3b	184
35be7ed7	185	AliMUONRecoParam* recoParam = 0x0;
6f6a3f3b	186
cc727848	187	AliMUONClusterStoreV2 clusterStore;
	188	AliMUONVClusterFinder* clusterFinder = 0x0;
	189
	190	AliMUONSimpleClusterServer* clusterServer = 0x0;
6f6a3f3b	191
cc727848	192	Int_t firstChamber(0);
cc727848	193	Int_t lastChamber(9);
6f6a3f3b	194
cc727848	195	if(whatToTest!=kOnlyTrigger){
	196	// Import TGeo geometry
	197	TString geoFilename = baseDir.Data();
	198	geoFilename.Append("/geometry.root");
	199	if ( ! AliGeomManager::GetGeometry() ) {
	200	AliGeomManager::LoadGeometry(geoFilename);
	201	if (! AliGeomManager::GetGeometry() ) {
6f6a3f3b	202	printf(">>> Error : getting geometry from file %s failed\n", geoFilename.Data());
6f6a3f3b	203	return;
cc727848	204	}
	205	}
	206	transformer = new AliMUONGeometryTransformer();
	207	// Load geometry data
	208	transformer->LoadGeometryData();
35be7ed7	209	// Load reconstruction parameters
	210	AliCDBPath path("MUON","Calib","RecoParam");
	211	AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
	212	if(entry) {
	213	recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
	214	entry->SetOwner(0);
	215	AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
	216	}
	217	if (!recoParam) {
	218	printf("Couldn't find RecoParam object in OCDB: create default one");
	219	recoParam = AliMUONRecoParam::GetLowFluxParam();
	220	}
	221	recoParam->Print("FULL");
cc727848	222	clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
0f211038	223	clusterServer = new AliMUONSimpleClusterServer(clusterFinder,*transformer);
cc727848	224	}
6f6a3f3b	225
cc727848	226	AliMUONDigitStoreV1 digitStore;
cc727848	227	AliMUONVDigit* mDigit;
6f6a3f3b	228
cc727848	229	Int_t clusterSize;
6f6a3f3b	230
cc727848	231	AliMUONDigitMaker digitMaker;
cc727848	232	TList digitsList[kNplanes];
6f6a3f3b	233
cc727848	234	Int_t nevents = RunLoader->GetNumberOfEvents();
cc727848	235	Int_t analysisFrac = nevents/10+1;
6f6a3f3b	236
cc727848	237	printf("\nNumber of events = %i\n\n",nevents);
6f6a3f3b	238
cc727848	239	for(Int_t ievent=0; ievent<nevents; ievent++){
	240	if(ievent%analysisFrac==0) printf("Analysing event = %i\n",ievent);
	241	RunLoader->GetEvent(ievent);
	242	if(whatToTest!=kOnlyTracker){
	243	digitStore.Clear();
	244	treeR = MUONLoader->TreeR();
	245	trigStore = AliMUONVTriggerStore::Create(*treeR);
6f6a3f3b	246
cc727848	247	if ( trigStore == 0x0 ) continue;
	248	trigStore->Clear();
	249	trigStore->Connect(*treeR);
	250	treeR->GetEvent(0);
6f6a3f3b	251
cc727848	252	TIter nextLocal(trigStore->CreateLocalIterator());
	253	while ( (locTrg = static_cast<AliMUONLocalTrigger*>( nextLocal() )) != NULL )
	254	{
6f6a3f3b	255	TArrayS xyPattern[2];
	256	locTrg->GetXPattern(xyPattern[0]);
	257	locTrg->GetYPattern(xyPattern[1]);
	258
	259	Int_t nBoard = locTrg->LoCircuit();
	260	digitMaker.TriggerDigits(nBoard, xyPattern, digitStore);
cc727848	261	}
6f6a3f3b	262
cc727848	263	TIter next(digitStore.CreateIterator());
6f6a3f3b	264
cc727848	265	while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
cc727848	266	{
6f6a3f3b	267	Int_t detElemId = mDigit->DetElemId();
	268	Int_t ch = detElemId/100 - 11;
	269	Int_t cathode = mDigit->Cathode();
	270	Int_t slat = detElemId%100;
	271	Int_t iplane = GetPlane(ch, cathode);
	272	trigStripMult[iplane]->Fill(slat);
	273	digitsList[iplane].Add(mDigit);
cc727848	274	} // loop on digits
6f6a3f3b	275
cc727848	276	for(Int_t iplane=0; iplane<kNplanes; iplane++){
6f6a3f3b	277	while(digitsList[iplane].GetEntries()){
	278	clusterSize=1;
	279	mDigit = (AliMUONVDigit*)digitsList[iplane].At(0);
	280	digitsList[iplane].Remove(mDigit);
	281	ClusterSize(digitsList[iplane],mDigit,clusterSize,1);
	282	//if(clusterSize>1) printf("Cluster size = %i\n\n",clusterSize);
	283	trigClusterSize[iplane]->Fill(clusterSize);
	284
	285	Int_t detElemId = mDigit->DetElemId();
	286	Int_t slat = detElemId%100;
	287	trigClusterMult[iplane]->Fill(slat);
	288	} // loop o sorted digits
cc727848	289	} // loop on planes
cc727848	290	} // trigger part
6f6a3f3b	291
cc727848	292	if(whatToTest!=kOnlyTrigger){
	293	clusterStore.Clear();
	294	treeD = MUONLoader->TreeD();
	295	digitStoreTrack = AliMUONVDigitStore::Create(*treeD);
6f6a3f3b	296
cc727848	297	if ( digitStoreTrack == 0x0 ) continue;
cc727848	298	digitStoreTrack->Clear();
6f6a3f3b	299	digitStoreTrackCut.Clear();
cc727848	300	digitStoreTrack->Connect(*treeD);
cc727848	301	treeD->GetEvent(0);
6f6a3f3b	302
	303	// Cut low charge channels (pedestal subtraction)
	304	TIter nextDigitTrack(digitStoreTrack->CreateIterator());
	305	while ( ( mDigit = static_cast<AliMUONVDigit*>(nextDigitTrack()) ) )
	306	{
	307	//printf("Digit class = %s",mDigit->ClassName());
	308	Float_t charge = mDigit->Charge();
	309	if(charge<adcCut) continue;
	310	digitStoreTrackCut.Add(*mDigit, AliMUONVDigitStore::kDeny);
	311	} // loop on digits
	312
0f211038	313	TIter nextDigitTrackCut(digitStoreTrackCut.CreateIterator());
7deb8eb0	314	clusterServer->UseDigits(nextDigitTrackCut);
6f6a3f3b	315
cc727848	316	for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
35be7ed7	317	clusterServer->Clusterize(ch, clusterStore, AliMpArea(),recoParam);
cc727848	318	}
6f6a3f3b	319
cc727848	320	TIter next(clusterStore.CreateIterator());
	321	while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) )
	322	{
6f6a3f3b	323	Float_t charge = cluster->GetCharge();
	324	if(charge==0.) continue;
	325	Int_t ch = cluster->GetChamberId();
	326	Int_t npads = cluster->GetNDigits();
	327	Int_t detElemId = cluster->GetDetElemId()%100;
	328	//printf("charge = %f pads = %i detElemId = %i\n",charge,npads,detElemId);
	329
	330	trackClusterSize[ch]->Fill(npads);
	331	trackClusterMult[ch]->Fill(charge);
	332	trackClusterYield[ch]->Fill(detElemId);
cc727848	333	} // loop on clusters
	334	} // tracker part
	335	} // loop on events
6f6a3f3b	336
cc727848	337	MUONLoader->UnloadRecPoints();
	338	MUONLoader->UnloadDigits();
	339	RunLoader->UnloadAll();
	340	delete RunLoader;
	341	TString outFileName = outDir.Data();
	342	outFileName.Append("/outTestRecPoints.root");
	343	TFile* outFile = new TFile(outFileName.Data(), "RECREATE");
	344	TDirectory* dir = 0x0;
	345	if(whatToTest!=kOnlyTracker){
	346	outFile->cd();
	347	dir = outFile->mkdir("Trigger");
	348	dir->cd();
	349	histoListTrig.Write();
	350	}
	351	if(whatToTest!=kOnlyTrigger){
	352	outFile->cd();
	353	dir = outFile->mkdir("Tracker");
	354	dir->cd();
	355	histoListTrack.Write();
	356	}
	357	outFile->Close();
	358	printf("\nSee results in %s\n",outFileName.Data());
	359	}
	360
	361	void ClusterSize(TList& list, AliMUONVDigit* refDigit, Int_t& clusterSize, Int_t depthLevel)
	362	{
	363	AliMUONVDigit* mDigit = 0x0;
	364	for(Int_t idigit=0; idigit<list.GetEntries(); idigit++){
	365	mDigit = (AliMUONVDigit*) list.At(idigit);
	366	if(mDigit->DetElemId() != refDigit->DetElemId()) continue;
	367	Int_t diffX = TMath::Abs(mDigit->PadX() - refDigit->PadX());
	368	Int_t diffY = TMath::Abs(mDigit->PadY() - refDigit->PadY());
	369	if(diffX>1) continue;
	370	if(diffY>1) continue;
	371	if(diffX*diffY != 0) continue;
	372	clusterSize++;
	373	list.Remove(mDigit);
	374	//printf("DetElemId = %i Level = %i Ref. (%2i,%2i) Found (%2i,%2i)\n",mDigit->DetElemId(),depthLevel,refDigit->PadX(),refDigit->PadY(),mDigit->PadX(),mDigit->PadY());
	375	ClusterSize(list, mDigit, clusterSize, depthLevel+1);
	376	Int_t val = idigit + depthLevel - clusterSize;
	377	//printf("Level = %i Current digit = %i\t",depthLevel,idigit);
	378	idigit = TMath::Max(-1,val);
	379	//printf("New digit = %i\n",idigit);
	380	}
	381	}