PWG0/trigger/AliGenTriggerMapSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AliGenTriggerMapSelector.h"
   4
   5 #include <TVector3.h>
   6 #include <TFile.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TNtuple.h>
  10 #include <TProfile.h>
  11 #include <TParticle.h>
  12 #include <TCanvas.h>
  13 #include <TTree.h>
  14 #include <TGeoManager.h>
  15 #include <TF1.h>
  16 #include <TGraph.h>
  17 #include <TLegend.h>
  18 #include <TLine.h>
  19
  20 #include <AliLog.h>
  21 #include <AliESD.h>
  22 #include <AliRunLoader.h>
  23 #include <AliStack.h>
  24 #include <AliGenEventHeader.h>
  25 #include <AliHeader.h>
  26
  27 #include <AliITSgeom.h>
  28 #include <AliITSLoader.h>
  29 #include <AliITSdigitSPD.h>
  30 #include <AliITSRecPoint.h>
  31
  32 #include "AliPWG0Helper.h"
  33
  34 //
  35 //
  36
  37 ClassImp(AliGenTriggerMapSelector)
  38
  39 AliGenTriggerMapSelector::AliGenTriggerMapSelector() :
  40   AliSelectorRL(),
  41   fChipsFired(0),
  42   fTracklets(0)
  43 {
  44   //
  45   // Constructor. Initialization of pointers
  46   //
  47 }
  48
  49 AliGenTriggerMapSelector::~AliGenTriggerMapSelector()
  50 {
  51   //
  52   // Destructor
  53   //
  54
  55   // histograms are in the output list and deleted when the output
  56   // list is deleted by the TSelector dtor
  57 }
  58
  59 void AliGenTriggerMapSelector::SlaveBegin(TTree *tree)
  60 {
  61   AliSelectorRL::SlaveBegin(tree);
  62
  63   fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5);
  64   fTracklets = new TNtuple("fTracklets", "", "vertex:chip1:chip2:dirty");
  65 }
  66
  67 void AliGenTriggerMapSelector::Init(TTree* tree)
  68 {
  69   // read the user objects
  70
  71   AliSelectorRL::Init(tree);
  72
  73   // enable only the needed branches
  74   if (tree)
  75   {
  76     tree->SetBranchStatus("*", 0);
  77     tree->SetBranchStatus("fTriggerMask", 1);
  78     tree->SetBranchStatus("fSPDVertex*", 1);
  79     tree->SetBranchStatus("fSPDMult*", 1);
  80   }
  81 }
  82
  83 Bool_t AliGenTriggerMapSelector::Process(Long64_t entry)
  84 {
  85   //
  86   // processing
  87   //
  88
  89   if (AliSelectorRL::Process(entry) == kFALSE)
  90     return kFALSE;
  91
  92   // Check prerequisites
  93   if (!fESD)
  94   {
  95     AliDebug(AliLog::kError, "ESD branch not available");
  96     return kFALSE;
  97   }
  98
  99   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);
 100
 101   if (!eventTriggered)
 102   {
 103     AliDebug(AliLog::kDebug, "Event not triggered. Skipping.");
 104     return kTRUE;
 105   }
 106
 107   AliHeader* header = GetHeader();
 108   if (!header)
 109   {
 110     AliDebug(AliLog::kError, "Header not available");
 111     return kFALSE;
 112   }
 113
 114   // get the MC vertex
 115   AliGenEventHeader* genHeader = header->GenEventHeader();
 116   TArrayF vtxMC(3);
 117   genHeader->PrimaryVertex(vtxMC);
 118   Float_t zVertex = vtxMC[2];
 119
 120   AliStack* stack = GetStack();
 121   if (!stack)
 122   {
 123     AliDebug(AliLog::kError, "Stack not available");
 124     return kFALSE;
 125   }
 126   Int_t nPrim  = stack->GetNprimary();
 127
 128   AliRunLoader* runLoader = GetRunLoader();
 129   if (!runLoader)
 130   {
 131     AliDebug(AliLog::kError, "runloader not available");
 132     return kFALSE;
 133   }
 134
 135   // TDirectory::TContext restores the current directory is restored when the scope ends.
 136   // This helps around ROOT bug #26025 and is good behaviour anyway
 137   TDirectory::TContext context(0);
 138   AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" );
 139   loader->LoadDigits("READ");
 140   TTree* treeD = loader->TreeD();
 141   if (!treeD)
 142   {
 143     AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS");
 144     return kFALSE;
 145   }
 146
 147   treeD->SetBranchStatus("*", 0);
 148   treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1);
 149   treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1);
 150
 151   TClonesArray* digits = 0;
 152   treeD->SetBranchAddress("ITSDigitsSPD", &digits);
 153   if (digits);
 154     digits->Clear();
 155
 156   const Int_t startSPD = 0; //geom->GetStartSPD();
 157   const Int_t lastSPD  = 239; //geom->GetLastSPD();
 158
 159   const AliMultiplicity* mult = fESD->GetMultiplicity();
 160   if (!mult)
 161   {
 162     AliDebug(AliLog::kError, "AliMultiplicity not available");
 163     return kFALSE;
 164   }
 165
 166   // loop over tracklets
 167   for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 168   {
 169     Float_t theta = mult->GetTheta(i);
 170     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
 171
 172     Int_t label = mult->GetLabel(i);
 173     if (label >= nPrim)
 174     {
 175       AliDebug(AliLog::kDebug, Form("Skipping particle %d because it is not a primary (label %d)", i, label));
 176       continue;
 177     }
 178     if (label < 0)
 179     {
 180       AliDebug(AliLog::kDebug, Form("Skipping particle %d because its label is negative (label %d)", i, label));
 181       continue;
 182     }
 183
 184     Int_t chip[2];
 185     Int_t nClusters[2];
 186
 187     for (Int_t j=0; j<2; ++j)
 188     {
 189       chip[j] = -1;
 190       nClusters[j] = 0;
 191     }
 192
 193     // now find clusters that belong to this label
 194     // loop over modules (ladders)
 195     for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++)
 196     {
 197       Int_t currentLayer = 0;
 198       if (moduleIndex >= 80)
 199         currentLayer = 1;
 200
 201       treeD->GetEvent(moduleIndex);
 202
 203       // get number of digits in this module
 204       Int_t ndigitsInModule = digits->GetEntriesFast();
 205
 206       // loop over digits in this module
 207       for (Int_t iDig=0; iDig<ndigitsInModule; iDig++)
 208       {
 209         AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig);
 210
 211         // check if digit belongs to this tracklet
 212         Bool_t belongs = kFALSE;
 213         for (Int_t iTrack = 0; iTrack < dp->GetNTracks(); ++iTrack)
 214         {
 215           //Printf("%d %d %d %d: %d =? %d", i, moduleIndex, iDig, iTrack, label, dp->GetTrack(iTrack));
 216           if (dp->GetTrack(iTrack) == label)
 217             belongs = kTRUE;
 218         }
 219
 220         if (!belongs)
 221           continue;
 222
 223         Int_t column = dp->GetCoord1();
 224         Int_t isChip = column / 32;
 225
 226         //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip);
 227
 228         fChipsFired->Fill(moduleIndex, isChip);
 229
 230         if (chip[currentLayer] == -1)
 231         {
 232           chip[currentLayer] = moduleIndex * 5 + isChip;
 233           //Printf("%d --> %d %d", label, moduleIndex, isChip);
 234           nClusters[currentLayer]++;
 235         }
 236         else
 237         {
 238           if (chip[currentLayer] != moduleIndex * 5 + isChip)
 239             nClusters[currentLayer]++;
 240         }
 241       }
 242     }
 243
 244     AliDebug(AliLog::kDebug, Form("Tracklet %d fired by chip %d and %d", i, chip[0], chip[1]));
 245
 246     Bool_t dirty = kFALSE;
 247     if (nClusters[0] == 0 || nClusters[1] == 0)
 248     {
 249       AliDebug(AliLog::kDebug, Form("Not enough clusters found for tracklet %d", i));
 250       continue;
 251     }
 252     if (nClusters[0] > 1 || nClusters[1] > 1)
 253     {
 254       AliDebug(AliLog::kDebug, Form("Too many clusters found for tracklet %d; L1: %d; L2: %d", i, nClusters[0], nClusters[1]));
 255       dirty = kTRUE;
 256     }
 257     fTracklets->Fill(zVertex, chip[0], chip[1], dirty);
 258   }
 259
 260   return kTRUE;
 261 }
 262
 263 Bool_t AliGenTriggerMapSelector::Notify()
 264 {
 265   AliRunLoader* runLoader = GetRunLoader();
 266
 267   if (runLoader)
 268   {
 269     AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" );
 270     if (loader)
 271     {
 272       loader->UnloadDigits();
 273       loader->UnloadRecPoints();
 274     }
 275   }
 276
 277   return AliSelectorRL::Notify();
 278 }
 279
 280 void AliGenTriggerMapSelector::SlaveTerminate()
 281 {
 282   // The SlaveTerminate() function is called after all entries or objects
 283   // have been processed. When running with PROOF SlaveTerminate() is called
 284   // on each slave server.
 285
 286   AliSelectorRL::SlaveTerminate();
 287
 288   // Add the histograms to the output on each slave server
 289   if (!fOutput)
 290   {
 291     AliDebug(AliLog::kError, "ERROR: Output list not initialized.");
 292     return;
 293   }
 294
 295   fOutput->Add(fChipsFired);
 296   fOutput->Add(fTracklets);
 297 }
 298
 299 void AliGenTriggerMapSelector::Terminate()
 300 {
 301   // The Terminate() function is the last function to be called during
 302   // a query. It always runs on the client, it can be used to present
 303   // the results graphically or save the results to file.
 304
 305   AliSelectorRL::Terminate();
 306
 307   fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired"));
 308   fTracklets = dynamic_cast<TNtuple*> (fOutput->FindObject("fTracklets"));
 309
 310   if (!fTracklets || !fChipsFired)
 311   {
 312     AliError("Histograms not available");
 313     return;
 314   }
 315
 316   WriteHistograms();
 317 }
 318
 319 void AliGenTriggerMapSelector::WriteHistograms(const char* filename)
 320 {
 321   TFile* file = TFile::Open(filename, "RECREATE");
 322
 323   fChipsFired->Write();
 324   fTracklets->Write();
 325
 326   file->Close();
 327 }
 328
 329 void AliGenTriggerMapSelector::ReadHistograms(const char* filename)
 330 {
 331   TFile* file = TFile::Open(filename);
 332
 333   if (!file)
 334     return;
 335
 336   fTracklets  = dynamic_cast<TNtuple*> (file->Get("fTracklets"));
 337   fChipsFired  = dynamic_cast<TH2F*> (file->Get("fChipsFired"));
 338 }
 339
 340 void AliGenTriggerMapSelector::GenerateTriggerMap(Bool_t clean)
 341 {
 342   Int_t nEntries = fTracklets->GetEntries();
 343   //nEntries = 10;
 344
 345   TH2F* hist = new TH2F("hist", ";layer1;layer2", 400, -0.5, 399.5, 800, 399.5, 1199.5);
 346
 347   Long64_t nSkipped = 0;
 348   for (Int_t i = 0; i < nEntries; ++i)
 349   {
 350     fTracklets->GetEntry(i);
 351     Float_t* vars = fTracklets->GetArgs();
 352
 353     if (clean && (Bool_t) vars[3] != kFALSE)
 354     {
 355       nSkipped++;
 356       continue;
 357     }
 358
 359     hist->Fill(vars[1], vars[2]);
 360   }
 361   if (clean)
 362     Printf("Skipped %lld dirty entries", nSkipped);
 363
 364   hist->Draw("COLZ");
 365
 366   for (Int_t i = 1; i <= hist->GetNbinsX(); ++i)
 367     for (Int_t j = 1; j <= hist->GetNbinsY(); ++j)
 368       if (hist->GetBinContent(i, j) > 0)
 369         Printf("(%d,%d) ", TMath::Nint(hist->GetXaxis()->GetBinCenter(i)), TMath::Nint(hist->GetYaxis()->GetBinCenter(j)));
 370 }