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
  42 // tracker
  43 #include "AliGeomManager.h"
  44
  45 // MUON
  46 #include "AliMpCDB.h"
  47 #include "AliMpDDL.h"
  48 #include "AliMpDDLStore.h"
  49 #include "AliMpConstants.h"
  50
  51 // trigger
  52 #include "AliMUONVTriggerStore.h"
  53 #include "AliMUONDigitStoreV1.h"
  54 #include "AliMUONDigitMaker.h"
  55 #include "AliMUONLocalTrigger.h"
  56
  57 // tracker
  58 #include "AliMUONClusterStoreV2.h"
  59 #include "AliMUONClusterFinderMLEM.h"
  60 #include "AliMUONPreClusterFinder.h"
  61 #include "AliMUONSimpleClusterServer.h"
  62 #include "AliMUONGeometryTransformer.h"
  63 #include "AliMUONVDigitStore.h"
  64 #include "AliMUONVCluster.h"
  65
  66 #endif
  67
  68 const Int_t kNtrigChambers = AliMpConstants::NofTriggerChambers();
  69 const Int_t kNcathodes = AliMpConstants::NofCathodes();
  70
  71 const Int_t kNtrackChambers = AliMpConstants::NofTrackingChambers();
  72
  73 enum {kOnlyTracker, kOnlyTrigger, kTrackTrig};
  74
  75 Int_t GetPlane(Int_t ch, Int_t cath){return kNcathodes * ch + cath;}
  76 void ClusterSize(TList&, AliMUONVDigit*, Int_t&, Int_t);
  77
  78 // Main Method
  79 void TestRecPoints(TString baseDir=".", TString outDir=".", Int_t whatToTest=kTrackTrig, Int_t runNumber=0, TString cdbStorage="local://$ALICE_ROOT")
  80 {
  81   const Int_t kNplanes = kNtrigChambers * kNcathodes;
  82   const Int_t kNslats = 18;
  83
  84   TList histoListTrig;
  85
  86   TH1F* trigStripMult[kNplanes];
  87   TH1F* trigClusterSize[kNplanes];
  88   TH1F* trigClusterMult[kNplanes];
  89   TString histoName, histoTitle;
  90
  91   if(whatToTest!=kOnlyTracker){
  92     for(Int_t ch=0; ch<kNtrigChambers; ch++){
  93       for(Int_t cath=0; cath<kNcathodes; cath++){
  94         Int_t iplane = GetPlane(ch, cath);
  95         histoName = "stripMult";
  96         histoName.Append(Form("Ch%iCath%i",ch,cath));
  97         histoTitle = Form("Strip multiplicity Ch %i Cathode %i", ch, cath);
  98         trigStripMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
  99         trigStripMult[iplane]->SetXTitle("slat");
 100         histoListTrig.Add(trigStripMult[iplane]);
 101
 102         histoName = "clusterSize";
 103         histoName.Append(Form("Ch%iCath%i",ch,cath));
 104         histoTitle = Form("Cluster size Ch %i Cathode %i", ch, cath);
 105         trigClusterSize[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), 10, 0.-0.5, (Float_t)10 - 0.5);
 106         trigClusterSize[iplane]->SetXTitle("cluster size");
 107         histoListTrig.Add(trigClusterSize[iplane]);
 108
 109         histoName = "clusterMult";
 110         histoName.Append(Form("Ch%iCath%i",ch,cath));
 111         histoTitle = Form("Cluster multiplicity Ch %i Cathode %i", ch, cath);
 112         trigClusterMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
 113         trigClusterMult[iplane]->SetXTitle("slat");
 114         histoListTrig.Add(trigClusterMult[iplane]);
 115       }
 116     }
 117   }
 118
 119   const Int_t kMaxDetElem = 26;
 120   TList histoListTrack;
 121   TH1F* trackClusterSize[kNtrackChambers];
 122   TH1F* trackClusterMult[kNtrackChambers];
 123   TH1F* trackClusterYield[kNtrackChambers];
 124
 125   if(whatToTest!=kOnlyTrigger){
 126     for(Int_t ch=0; ch<kNtrackChambers; ch++){
 127       histoName = "clusterSize";
 128       histoName.Append(Form("Ch%i",ch));
 129       histoTitle = Form("Cluster size Ch %i", ch);
 130       trackClusterSize[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
 131       trackClusterSize[ch]->SetXTitle("cluster size");
 132       histoListTrack.Add(trackClusterSize[ch]);
 133
 134       histoName = "clusterMult";
 135       histoName.Append(Form("Ch%i",ch));
 136       histoTitle = Form("Cluster multiplicity vs. Charge Ch %i", ch);
 137       trackClusterMult[ch] = new TH1F(histoName.Data(), histoTitle.Data(), 200, 0, 15000.);
 138       trackClusterMult[ch]->SetXTitle("Charge (pC)");
 139       histoListTrack.Add(trackClusterMult[ch]);
 140
 141       histoName = "clusterYield";
 142       histoName.Append(Form("Ch%i",ch));
 143       histoTitle = Form("Cluster yield vs. DetElem Ch %i", ch);
 144       trackClusterYield[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
 145       trackClusterYield[ch]->SetXTitle("Detector element #");
 146       histoListTrack.Add(trackClusterYield[ch]);
 147     }
 148   }
 149
 150
 151   // Creating Run Loader and opening RecPoints
 152   TString filename = baseDir.Data();
 153   filename.Append("/galice.root");
 154   AliRunLoader * RunLoader = AliRunLoader::Open(filename.Data(),"MUONLoader","UPDATE");
 155   if (RunLoader ==0x0) {
 156     printf(">>> Error : Error Opening %s file \n",filename.Data());
 157     return;
 158   }
 159
 160   // Loading MUON subsystem
 161   AliLoader* MUONLoader = RunLoader->GetDetectorLoader("MUON");
 162   if(whatToTest!=kOnlyTracker)  MUONLoader->LoadRecPoints("READ");
 163   if(whatToTest!=kOnlyTrigger)  MUONLoader->LoadDigits("READ");
 164
 165   TTree *treeR = 0x0, *treeD = 0x0;
 166   AliMUONVTriggerStore* trigStore = 0x0;
 167   AliMUONLocalTrigger* locTrg = 0x0;
 168
 169   AliMUONVDigitStore* digitStoreTrack = 0x0;
 170   AliMUONVCluster* cluster = 0x0;
 171
 172   // Load segmentation
 173   AliCDBManager::Instance()->SetDefaultStorage(cdbStorage.Data());
 174   AliCDBManager::Instance()->SetRun(runNumber);
 175
 176   AliMpCDB::LoadDDLStore();
 177
 178   AliMUONGeometryTransformer* transformer = 0x0;
 179
 180
 181   AliMUONClusterStoreV2 clusterStore;
 182   AliMUONVClusterFinder* clusterFinder = 0x0;
 183
 184   AliMUONSimpleClusterServer* clusterServer = 0x0;
 185
 186   Int_t firstChamber(0);
 187   Int_t lastChamber(9);
 188
 189   if(whatToTest!=kOnlyTrigger){
 190     // Import TGeo geometry
 191     TString geoFilename = baseDir.Data();
 192     geoFilename.Append("/geometry.root");
 193     if ( ! AliGeomManager::GetGeometry() ) {
 194       AliGeomManager::LoadGeometry(geoFilename);
 195       if (! AliGeomManager::GetGeometry() ) {
 196         printf(">>> Error : getting geometry from file %s failed\n", geoFilename.Data());
 197         return;
 198       }
 199     }
 200     transformer = new AliMUONGeometryTransformer();
 201     // Load geometry data
 202     transformer->LoadGeometryData();
 203     clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
 204     clusterServer = new AliMUONSimpleClusterServer(*clusterFinder,*transformer);
 205   }
 206
 207   AliMUONDigitStoreV1 digitStore;
 208   AliMUONVDigit* mDigit;
 209
 210   Int_t clusterSize;
 211
 212   AliMUONDigitMaker digitMaker;
 213   TList digitsList[kNplanes];
 214
 215   Int_t nevents = RunLoader->GetNumberOfEvents();
 216   Int_t analysisFrac = nevents/10+1;
 217
 218   printf("\nNumber of events = %i\n\n",nevents);
 219
 220   for(Int_t ievent=0; ievent<nevents; ievent++){
 221     if(ievent%analysisFrac==0) printf("Analysing event = %i\n",ievent);
 222     RunLoader->GetEvent(ievent);
 223     if(whatToTest!=kOnlyTracker){
 224       digitStore.Clear();
 225       treeR = MUONLoader->TreeR();
 226       trigStore = AliMUONVTriggerStore::Create(*treeR);
 227
 228       if ( trigStore == 0x0 ) continue;
 229       trigStore->Clear();
 230       trigStore->Connect(*treeR);
 231       treeR->GetEvent(0);
 232
 233       TIter nextLocal(trigStore->CreateLocalIterator());
 234       while ( (locTrg = static_cast<AliMUONLocalTrigger*>( nextLocal() )) != NULL )
 235       {
 236         TArrayS xyPattern[2];
 237         locTrg->GetXPattern(xyPattern[0]);
 238         locTrg->GetYPattern(xyPattern[1]);
 239
 240         Int_t nBoard = locTrg->LoCircuit();
 241         digitMaker.TriggerDigits(nBoard, xyPattern, digitStore);
 242       }
 243
 244       TIter next(digitStore.CreateIterator());
 245
 246       while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
 247       {
 248         Int_t detElemId = mDigit->DetElemId();
 249         Int_t ch = detElemId/100 - 11;
 250         Int_t cathode = mDigit->Cathode();
 251         Int_t slat = detElemId%100;
 252         Int_t iplane = GetPlane(ch, cathode);
 253         trigStripMult[iplane]->Fill(slat);
 254         digitsList[iplane].Add(mDigit);
 255       } // loop on digits
 256
 257       for(Int_t iplane=0; iplane<kNplanes; iplane++){
 258         while(digitsList[iplane].GetEntries()){
 259           clusterSize=1;
 260           mDigit = (AliMUONVDigit*)digitsList[iplane].At(0);
 261           digitsList[iplane].Remove(mDigit);
 262           ClusterSize(digitsList[iplane],mDigit,clusterSize,1);
 263           //if(clusterSize>1) printf("Cluster size = %i\n\n",clusterSize);
 264           trigClusterSize[iplane]->Fill(clusterSize);
 265
 266           Int_t detElemId = mDigit->DetElemId();
 267           Int_t slat = detElemId%100;
 268           trigClusterMult[iplane]->Fill(slat);
 269         } // loop o sorted digits
 270       } // loop on planes
 271     } // trigger part
 272
 273     if(whatToTest!=kOnlyTrigger){
 274       clusterStore.Clear();
 275       treeD = MUONLoader->TreeD();
 276       digitStoreTrack = AliMUONVDigitStore::Create(*treeD);
 277
 278       if ( digitStoreTrack == 0x0 ) continue;
 279       digitStoreTrack->Clear();
 280       digitStoreTrack->Connect(*treeD);
 281       treeD->GetEvent(0);
 282       clusterServer->UseDigitStore(*digitStoreTrack);
 283
 284       for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
 285         clusterServer->Clusterize(ch, clusterStore, AliMpArea());
 286       }
 287
 288       TIter next(clusterStore.CreateIterator());
 289       while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) )
 290       {
 291         Float_t charge = cluster->GetCharge();
 292         if(charge==0.) continue;
 293         Int_t ch = cluster->GetChamberId();
 294         Int_t npads = cluster->GetNDigits();
 295         Int_t detElemId = cluster->GetDetElemId()%100;
 296         //printf("charge = %f   pads = %i   detElemId = %i\n",charge,npads,detElemId);
 297
 298         trackClusterSize[ch]->Fill(npads);
 299         trackClusterMult[ch]->Fill(charge);
 300         trackClusterYield[ch]->Fill(detElemId);
 301       } // loop on clusters
 302     } // tracker part
 303   } // loop on events
 304
 305   MUONLoader->UnloadRecPoints();
 306   MUONLoader->UnloadDigits();
 307   RunLoader->UnloadAll();
 308   delete RunLoader;
 309   TString outFileName = outDir.Data();
 310   outFileName.Append("/outTestRecPoints.root");
 311   TFile* outFile = new TFile(outFileName.Data(), "RECREATE");
 312   TDirectory* dir = 0x0;
 313   if(whatToTest!=kOnlyTracker){
 314     outFile->cd();
 315     dir = outFile->mkdir("Trigger");
 316     dir->cd();
 317     histoListTrig.Write();
 318   }
 319   if(whatToTest!=kOnlyTrigger){
 320     outFile->cd();
 321     dir = outFile->mkdir("Tracker");
 322     dir->cd();
 323     histoListTrack.Write();
 324   }
 325   outFile->Close();
 326   printf("\nSee results in %s\n",outFileName.Data());
 327 }
 328
 329 void ClusterSize(TList& list, AliMUONVDigit* refDigit, Int_t& clusterSize, Int_t depthLevel)
 330 {
 331   AliMUONVDigit* mDigit = 0x0;
 332   for(Int_t idigit=0; idigit<list.GetEntries(); idigit++){
 333     mDigit = (AliMUONVDigit*) list.At(idigit);
 334     if(mDigit->DetElemId() != refDigit->DetElemId()) continue;
 335     Int_t diffX = TMath::Abs(mDigit->PadX() - refDigit->PadX());
 336     Int_t diffY = TMath::Abs(mDigit->PadY() - refDigit->PadY());
 337     if(diffX>1) continue;
 338     if(diffY>1) continue;
 339     if(diffX*diffY != 0) continue;
 340     clusterSize++;
 341     list.Remove(mDigit);
 342     //printf("DetElemId = %i   Level = %i   Ref. (%2i,%2i)    Found (%2i,%2i)\n",mDigit->DetElemId(),depthLevel,refDigit->PadX(),refDigit->PadY(),mDigit->PadX(),mDigit->PadY());
 343     ClusterSize(list, mDigit, clusterSize, depthLevel+1);
 344     Int_t val = idigit + depthLevel - clusterSize;
 345     //printf("Level = %i   Current digit = %i\t",depthLevel,idigit);
 346     idigit = TMath::Max(-1,val);
 347     //printf("New digit = %i\n",idigit);
 348   }
 349 }