MUON/TestRecPoints.C

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $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"
  41 #include "AliCDBEntry.h"
  42 #include "AliCDBPath.h"
  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"
  66 #include "AliMUONDigitStoreV2R.h"
  67 #include "AliMUONVCluster.h"
  68 #include "AliMUONRecoParam.h"
  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
  83 void TestRecPoints(TString baseDir=".", TString outDir=".", Float_t adcCut = 10., Int_t whatToTest=kTrackTrig, Int_t runNumber=0, TString cdbStorage="local://$ALICE_ROOT/OCDB")
  84 {
  85   const Int_t kNplanes = kNtrigChambers * kNcathodes;
  86   const Int_t kNslats = 18;
  87
  88   TList histoListTrig;
  89
  90   TH1F* trigStripMult[kNplanes];
  91   TH1F* trigClusterSize[kNplanes];
  92   TH1F* trigClusterMult[kNplanes];
  93   TString histoName, histoTitle;
  94
  95   if(whatToTest!=kOnlyTracker){
  96     for(Int_t ch=0; ch<kNtrigChambers; ch++){
  97       for(Int_t cath=0; cath<kNcathodes; cath++){
  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]);
 119       }
 120     }
 121   }
 122
 123   const Int_t kMaxDetElem = 26;
 124   TList histoListTrack;
 125   TH1F* trackClusterSize[kNtrackChambers];
 126   TH1F* trackClusterMult[kNtrackChambers];
 127   TH1F* trackClusterYield[kNtrackChambers];
 128
 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]);
 137
 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.);
 142       trackClusterMult[ch]->SetXTitle("Charge (ADC)");
 143       histoListTrack.Add(trackClusterMult[ch]);
 144
 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
 154
 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   }
 163
 164   // Loading MUON subsystem
 165   AliLoader* MUONLoader = RunLoader->GetDetectorLoader("MUON");
 166   if(whatToTest!=kOnlyTracker)  MUONLoader->LoadRecPoints("READ");
 167   if(whatToTest!=kOnlyTrigger)  MUONLoader->LoadDigits("READ");
 168
 169   TTree *treeR = 0x0, *treeD = 0x0;
 170   AliMUONVTriggerStore* trigStore = 0x0;
 171   AliMUONLocalTrigger* locTrg = 0x0;
 172
 173   AliMUONVDigitStore* digitStoreTrack = 0x0;
 174   AliMUONDigitStoreV2R digitStoreTrackCut;
 175   AliMUONVCluster* cluster = 0x0;
 176
 177   // Load segmentation
 178   AliCDBManager::Instance()->SetDefaultStorage(cdbStorage.Data());
 179   AliCDBManager::Instance()->SetRun(runNumber);
 180
 181   AliMpCDB::LoadDDLStore();
 182
 183   AliMUONGeometryTransformer* transformer = 0x0;
 184
 185   AliMUONRecoParam* recoParam = 0x0;
 186
 187   AliMUONClusterStoreV2 clusterStore;
 188   AliMUONVClusterFinder* clusterFinder = 0x0;
 189
 190   AliMUONSimpleClusterServer* clusterServer = 0x0;
 191
 192   Int_t firstChamber(0);
 193   Int_t lastChamber(9);
 194
 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() ) {
 202         printf(">>> Error : getting geometry from file %s failed\n", geoFilename.Data());
 203         return;
 204       }
 205     }
 206     transformer = new AliMUONGeometryTransformer();
 207     // Load geometry data
 208     transformer->LoadGeometryData();
 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");
 222     clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
 223     clusterServer = new AliMUONSimpleClusterServer(clusterFinder,*transformer);
 224   }
 225
 226   AliMUONDigitStoreV1 digitStore;
 227   AliMUONVDigit* mDigit;
 228
 229   Int_t clusterSize;
 230
 231   AliMUONDigitMaker digitMaker;
 232   TList digitsList[kNplanes];
 233
 234   Int_t nevents = RunLoader->GetNumberOfEvents();
 235   Int_t analysisFrac = nevents/10+1;
 236
 237   printf("\nNumber of events = %i\n\n",nevents);
 238
 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);
 246
 247       if ( trigStore == 0x0 ) continue;
 248       trigStore->Clear();
 249       trigStore->Connect(*treeR);
 250       treeR->GetEvent(0);
 251
 252       TIter nextLocal(trigStore->CreateLocalIterator());
 253       while ( (locTrg = static_cast<AliMUONLocalTrigger*>( nextLocal() )) != NULL )
 254       {
 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);
 261       }
 262
 263       TIter next(digitStore.CreateIterator());
 264
 265       while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
 266       {
 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);
 274       } // loop on digits
 275
 276       for(Int_t iplane=0; iplane<kNplanes; iplane++){
 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
 289       } // loop on planes
 290     } // trigger part
 291
 292     if(whatToTest!=kOnlyTrigger){
 293       clusterStore.Clear();
 294       treeD = MUONLoader->TreeD();
 295       digitStoreTrack = AliMUONVDigitStore::Create(*treeD);
 296
 297       if ( digitStoreTrack == 0x0 ) continue;
 298       digitStoreTrack->Clear();
 299       digitStoreTrackCut.Clear();
 300       digitStoreTrack->Connect(*treeD);
 301       treeD->GetEvent(0);
 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
 313       TIter nextDigitTrackCut(digitStoreTrackCut.CreateIterator());
 314       clusterServer->UseDigits(nextDigitTrackCut,&digitStoreTrackCut);
 315
 316       for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
 317         clusterServer->Clusterize(ch, clusterStore, AliMpArea(),recoParam);
 318       }
 319
 320       TIter next(clusterStore.CreateIterator());
 321       while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) )
 322       {
 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);
 333       } // loop on clusters
 334     } // tracker part
 335   } // loop on events
 336
 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 }