PWG0/TPC/AliROCESDAnalysisSelector.cxx

   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$ */
  17
  18 //
  19 //  This class analyses TPC cosmics data from the ESD and the ESDfriend
  20 //
  21 //  Authors: Jan.Fiete.Grosse-Oetringhaus@cern.ch, Claus.Jorgensen@cern.ch
  22 //
  23
  24 #include "AliROCESDAnalysisSelector.h"
  25
  26 #include <AliLog.h>
  27 #include <AliESD.h>
  28 #include <AliESDfriend.h>
  29 #include <AliTPCclusterMI.h>
  30 #include <AliTPCseed.h>
  31
  32 #include <TFile.h>
  33 #include <TMath.h>
  34 #include <TTree.h>
  35 #include <TCanvas.h>
  36 #include <TSystem.h>
  37 #include <TObjArray.h>
  38 #include <TTimeStamp.h>
  39
  40 #include "TPC/AliTPCClusterHistograms.h"
  41
  42 extern TSystem* gSystem;
  43
  44 ClassImp(AliROCESDAnalysisSelector)
  45
  46 AliROCESDAnalysisSelector::AliROCESDAnalysisSelector() :
  47   AliSelector(),
  48   fESDfriend(0),
  49   fObjectsToSave(0)
  50 {
  51   //
  52   // Constructor. Initialization of pointers
  53   //
  54   fMinNumberOfRowsIsTrack = 5;
  55
  56   fObjectsToSave = new TObjArray();
  57
  58   for (Int_t i=0; i<kTPCHists; i++)
  59     fClusterHistograms[i] = 0;
  60 }
  61
  62 AliROCESDAnalysisSelector::~AliROCESDAnalysisSelector()
  63 {
  64   //
  65   // Destructor
  66   //
  67 }
  68
  69 void AliROCESDAnalysisSelector::SlaveBegin(TTree* tree)
  70 {
  71   //
  72
  73   AliSelector::SlaveBegin(tree);
  74 }
  75
  76 void AliROCESDAnalysisSelector::Init(TTree *tree)
  77 {
  78   // The Init() function is called when the selector needs to initialize
  79   // a new tree or chain. Typically here the branch addresses of the tree
  80   // will be set. It is normaly not necessary to make changes to the
  81   // generated code, but the routine can be extended by the user if needed.
  82   // Init() will be called many times when running with PROOF.
  83
  84   AliSelector::Init(tree);
  85
  86   printf("Init called %p\n", (void*) fESDfriend);
  87
  88   // Set branch address
  89   if (tree)
  90   {
  91     tree->SetBranchAddress("ESDfriend", &fESDfriend);
  92
  93     tree->SetBranchStatus("*", 0);
  94     tree->SetBranchStatus("fTracks.*", 1);
  95     tree->SetBranchStatus("fTimeStamp", 1);
  96     //tree->SetBranchStatus("fTracks.fCalibContainer", 0);
  97   }
  98
  99   if (fESDfriend != 0)
 100     AliDebug(AliLog::kInfo, "INFO: Found ESDfriend branch in chain.");
 101 }
 102
 103 Bool_t AliROCESDAnalysisSelector::Process(Long64_t entry)
 104 {
 105   //
 106   // Implement your analysis here. Do not forget to call the parent class Process by
 107   // if (AliSelector::Process(entry) == kFALSE)
 108   //   return kFALSE;
 109   //
 110
 111   if (AliSelector::Process(entry) == kFALSE)
 112     return kFALSE;
 113
 114   // Check prerequisites
 115   if (!fESD)
 116   {
 117     AliDebug(AliLog::kError, "ESD branch not available");
 118     return kFALSE;
 119   }
 120
 121   // Check prerequisites
 122   if (!fESDfriend)
 123   {
 124     AliDebug(AliLog::kError, "ESDfriend branch not available");
 125     return kFALSE;
 126   }
 127
 128   fESD->SetESDfriend(fESDfriend);
 129
 130   Int_t flag = ProcessEvent(entry, kFALSE);
 131   if (flag == 1)
 132     ProcessEvent(entry, kTRUE);
 133
 134   // TODO This should not be needed, the TTree::GetEntry() should take care of this, maybe because it has a reference member, to be analyzed
 135   // if the ESDfriend is not deleted we get a major memory leak
 136   // here the esdfriend seems to be also deleted, very weird behaviour....
 137
 138   delete fESD;
 139   fESD = 0;
 140
 141   //delete fESDfriend;
 142   //fESDfriend = 0;
 143
 144
 145   return kTRUE;
 146 }
 147
 148 Int_t AliROCESDAnalysisSelector::ProcessEvent(Long64_t entry, Bool_t detailedHistogram)
 149 {
 150   //
 151   // Looping over tracks and clusters in event and filling histograms
 152   //
 153   // - if detailedHistogram = kTRUE special histograms are saved (in fObjectsToSave)
 154   // - the method returns
 155   //   1 : if a "flash" is detected something special in this event
 156   //
 157
 158   // save maximum 50 objects
 159   if (detailedHistogram)
 160     if (fObjectsToSave->GetEntries() > 50)
 161       return 0;
 162
 163   // for saving single events
 164   AliTPCClusterHistograms* clusterHistograms[kTPCSectors];
 165   for (Int_t i=0; i<kTPCSectors; i++)
 166     clusterHistograms[i] = 0;
 167
 168   Bool_t intToReturn = 0;
 169
 170   Int_t nTracks = fESD->GetNumberOfTracks();
 171
 172   Int_t nSkippedSeeds = 0;
 173   Int_t nSkippedTracks = 0;
 174
 175   // for "flash" detection
 176   Int_t nClusters = 0;
 177   Float_t clusterQtotSumVsTime[250];
 178   for (Int_t z=0; z<250; z++)
 179     clusterQtotSumVsTime[z] = 0;
 180
 181   // loop over esd tracks
 182   for (Int_t t=0; t<nTracks; t++)
 183   {
 184     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (fESD->GetTrack(t));
 185     if (!esdTrack)
 186     {
 187       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", t));
 188       continue;
 189     }
 190
 191     AliESDfriendTrack* friendtrack = const_cast<AliESDfriendTrack*> (dynamic_cast<const AliESDfriendTrack*> (esdTrack->GetFriendTrack()));
 192     if (!friendtrack)
 193     {
 194       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve friend of track %d.", t));
 195       continue;
 196     }
 197
 198     const AliTPCseed* seed = dynamic_cast<const AliTPCseed*> (friendtrack->GetCalibObject(0));
 199     if (!seed)
 200     {
 201       AliDebug(AliLog::kDebug, Form("ERROR: Could not retrieve seed of track %d.", t));
 202       nSkippedSeeds++;
 203       continue;
 204     }
 205
 206     if (!AcceptTrack(seed, fMinNumberOfRowsIsTrack))
 207     {
 208       AliDebug(AliLog::kDebug, Form("INFO: Rejected track %d.", t));
 209       nSkippedTracks++;
 210       continue;
 211     }
 212
 213     for (Int_t clusterID = 0; clusterID < 160; clusterID++)
 214     {
 215       AliTPCclusterMI* cluster = seed->GetClusterPointer(clusterID);
 216       if (!cluster)
 217         continue;
 218
 219       Int_t detector = cluster->GetDetector();
 220
 221       if (detector < 0 || detector >= kTPCSectors)
 222       {
 223         AliDebug(AliLog::kDebug, Form("We found a cluster from invalid sector %d", detector));
 224         continue;
 225       }
 226
 227       if (!detailedHistogram) {
 228         // TODO: find a clever way to handle the time
 229         Int_t time = 0;
 230
 231         if (fESD->GetTimeStamp()>1160000000)
 232           time = fESD->GetTimeStamp();
 233
 234         if (!fClusterHistograms[detector])
 235           fClusterHistograms[detector] = new AliTPCClusterHistograms(detector,"",time,time+5*60*60);
 236
 237         if (!fClusterHistograms[detector+kTPCSectors])
 238           fClusterHistograms[detector+kTPCSectors] = new AliTPCClusterHistograms(detector,"",time,time+5*60*60, kTRUE);
 239
 240         fClusterHistograms[detector]->FillCluster(cluster, time);
 241         fClusterHistograms[detector+kTPCSectors]->FillCluster(cluster, time);
 242
 243         Int_t z = Int_t(cluster->GetZ());
 244         if (z>=0 && z<250) {
 245           nClusters++;
 246           clusterQtotSumVsTime[z] += cluster->GetQ();
 247         }
 248       } // end of if !detailedHistograms
 249       else {
 250         // if we need the detailed histograms for this event
 251         if (!clusterHistograms[detector])
 252           clusterHistograms[detector] = new AliTPCClusterHistograms(detector, Form("flash_entry%d", entry));
 253
 254         clusterHistograms[detector]->FillCluster(cluster);
 255       }
 256     }
 257
 258     for (Int_t i=0; i<kTPCHists; i++)
 259       if (fClusterHistograms[i])
 260         fClusterHistograms[i]->FillTrack(seed);
 261
 262   }
 263
 264   // check if there's a very large q deposit ("flash")
 265   if (!detailedHistogram) {
 266     for (Int_t z=0; z<250; z++) {
 267       if (clusterQtotSumVsTime[z] > 150000) {
 268         printf(Form("  \n   -> Entry %lld sum of clusters at time %d is %f, ESD timestamp: %s (%d) \n \n", entry, z, clusterQtotSumVsTime[z], TTimeStamp(fESD->GetTimeStamp()).AsString(), fESD->GetTimeStamp()));
 269         intToReturn = 1;
 270       }
 271     }
 272   }
 273   else {
 274     for (Int_t i=0; i< kTPCSectors; i++) {
 275       if (clusterHistograms[i]) {
 276         fObjectsToSave->Add(clusterHistograms[i]);
 277       }
 278     }
 279   }
 280
 281 //   if (nSkippedSeeds > 0)
 282 //     printf("WARNING: The seed was not found for %d out of %d tracks.\n", nSkippedSeeds, nTracks);
 283 //   if (nSkippedTracks > 0)
 284 //     printf("INFO: Rejected %d out of %d tracks.\n", nSkippedTracks, nTracks);
 285
 286   return intToReturn;
 287 }
 288
 289 Bool_t AliROCESDAnalysisSelector::AcceptTrack(const AliTPCseed* track, Int_t minRowsIncluded) {
 290   //
 291   // check if the track should be accepted.
 292   //
 293   const Int_t   kMinClusters = 5;
 294   const Float_t kMinRatio    = 0.75;
 295   const Float_t kMax1pt      = 0.5;
 296
 297   Int_t  nRowsUsedByTracks = 0;
 298   Bool_t rowIncluded[96];
 299
 300   Float_t totalQtot = 0;
 301   Int_t   nClusters = 0;
 302
 303   for(Int_t r=0; r<96; r++)
 304     rowIncluded[r] = kFALSE;
 305
 306   for (Int_t clusterID = 0; clusterID < 160; clusterID++) {
 307     AliTPCclusterMI* cluster = track->GetClusterPointer(clusterID);
 308
 309     if (!cluster)
 310       continue;
 311
 312     Float_t qTot =   cluster->GetQ();
 313
 314     nClusters++;
 315     totalQtot += qTot;
 316
 317     if (!rowIncluded[cluster->GetRow()]) {
 318       nRowsUsedByTracks++;
 319       rowIncluded[cluster->GetRow()] = kTRUE;
 320     }
 321   }
 322
 323   Float_t meanQtot = totalQtot/nClusters;
 324
 325   if (meanQtot<70)
 326     return kFALSE;
 327
 328   if (nRowsUsedByTracks < minRowsIncluded)
 329     return kFALSE;
 330
 331   //  printf(Form("    TRACK: n clusters = %d,  n pad rows = %d \n",track->GetNumberOfClusters(), nRowsUsedByTracks));
 332
 333   if (track->GetNumberOfClusters()<kMinClusters) return kFALSE;
 334   Float_t ratio = track->GetNumberOfClusters()/(track->GetNFoundable()+1.);
 335   if (ratio<kMinRatio) return kFALSE;
 336   Float_t mpt = track->Get1Pt();
 337   if (TMath::Abs(mpt)>kMax1pt) return kFALSE;
 338
 339   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
 340   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
 341   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
 342
 343   return kTRUE;
 344 }
 345
 346 void AliROCESDAnalysisSelector::SlaveTerminate()
 347 {
 348   //
 349
 350   if (fOutput)
 351   {
 352     for (Int_t i=0; i<kTPCHists; i++)
 353       if (fClusterHistograms[i])
 354         fOutput->Add(fClusterHistograms[i]);
 355   }
 356 }
 357
 358 void AliROCESDAnalysisSelector::Terminate()
 359 {
 360   //
 361   // read the objects from the output list and write them to a file
 362   // the filename is modified by the object comment passed in the tree info or input list
 363   //
 364
 365   if (fOutput)
 366   {
 367     fOutput->Print();
 368
 369     for (Int_t i=0; i<kTPCSectors; i++)
 370       fClusterHistograms[i] = dynamic_cast<AliTPCClusterHistograms*> (fOutput->FindObject(AliTPCClusterHistograms::FormDetectorName(i, kFALSE)));
 371     for (Int_t i=0; i<kTPCSectors; i++)
 372       fClusterHistograms[kTPCSectors+i] = dynamic_cast<AliTPCClusterHistograms*> (fOutput->FindObject(AliTPCClusterHistograms::FormDetectorName(i, kTRUE)));
 373   }
 374
 375   TNamed* comment = 0;
 376   if (fTree && fTree->GetUserInfo())
 377     comment = dynamic_cast<TNamed*>(fTree->GetUserInfo()->FindObject("comment"));
 378   if (!comment && fInput)
 379     comment = dynamic_cast<TNamed*>(fInput->FindObject("comment"));
 380
 381   if (comment)
 382   {
 383     AliDebug(AliLog::kInfo, Form("INFO: Found comment in input list: %s", comment->GetTitle()));
 384   }
 385   else
 386     return;
 387
 388   TFile* file = TFile::Open(Form("rocESD_%s.root",comment->GetTitle()), "RECREATE");
 389
 390   for (Int_t i=0; i<kTPCHists; i++)
 391     if (fClusterHistograms[i]) {
 392       fClusterHistograms[i]->SaveHistograms();
 393       TCanvas* c = fClusterHistograms[i]->DrawHistograms();
 394                         TString dir;
 395                         dir.Form("WWW/%s/%s", comment->GetTitle(), c->GetName());
 396                         gSystem->mkdir(dir, kTRUE);
 397       c->SaveAs(Form("%s/plots_%s_%s.eps",dir.Data(),comment->GetTitle(),c->GetName()));
 398       c->SaveAs(Form("%s/plots_%s_%s.gif",dir.Data(),comment->GetTitle(),c->GetName()));
 399
 400       c->Close();
 401       delete c;
 402     }
 403
 404   gDirectory->mkdir("saved_objects");
 405   gDirectory->cd("saved_objects");
 406
 407   for (Int_t i=0; i<fObjectsToSave->GetSize(); i++) {
 408     if (fObjectsToSave->At(i)) {
 409       AliTPCClusterHistograms* clusterHistograms = dynamic_cast<AliTPCClusterHistograms*> (fObjectsToSave->At(i));
 410       if (clusterHistograms)
 411         clusterHistograms->SaveHistograms();
 412       else
 413         fObjectsToSave->At(i)->Write();
 414     }
 415   }
 416
 417   gDirectory->cd("../");
 418
 419
 420   file->Close();
 421 }