]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/TPC/AliROCESDAnalysisSelector.cxx
Typo corrected.
[u/mrichter/AliRoot.git] / PWG0 / TPC / AliROCESDAnalysisSelector.cxx
index af3ffec475d775c00329ff044b6a2c1f531f1141..6328d53920f107ab72bdb8cd68576eb922fb00e4 100644 (file)
 
 /* $Id$ */
 
-// The ESD is available as member fESD
+// 
+//  This class analyses TPC cosmics data from the ESD and the ESDfriend
 //
-// The Process function is nearly empty. Implement your analysis there and look at the other listed below functions you
-// might need.
+//  Authors: Jan.Fiete.Grosse-Oetringhaus@cern.ch, Claus.Jorgensen@cern.ch
 //
-// The following methods can be overrriden. Please do not forgot to call the base class function.
-//
-//    Begin():        called everytime a loop on the tree starts,
-//                    a convenient place to create your histograms.
-//    SlaveBegin():   called after Begin(), when on PROOF called only on the
-//                    slave servers.
-//    Init():         called for each new tree. Enable/Disable branches here.
-//    Process():      called for each event, in this function you decide what
-//                    to read and fill your histograms.
-//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
-//                    called only on the slave servers.
-//    Terminate():    called at the end of the loop on the tree,
-//                    a convenient place to draw/fit your histograms.
-//
-//  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
 
 #include "AliROCESDAnalysisSelector.h"
 
 #include <../TPC/AliTPCseed.h>
 
 #include <TFile.h>
+#include <TMath.h>
 #include <TTree.h>
 #include <TCanvas.h>
+#include <TSystem.h>
+#include <TObjArray.h>
+#include <TTimeStamp.h>
 
 #include "TPC/AliTPCClusterHistograms.h"
 
+extern TSystem* gSystem;
+
 ClassImp(AliROCESDAnalysisSelector)
 
 AliROCESDAnalysisSelector::AliROCESDAnalysisSelector() :
   AliSelector(),
   fESDfriend(0),
-  fClusterHistograms(0)
+  fObjectsToSave(0)
 {
   //
   // Constructor. Initialization of pointers
   //
+  fMinNumberOfRowsIsTrack = 5;
+
+  fObjectsToSave = new TObjArray();
+  
+  for (Int_t i=0; i<kTPCHists; i++)
+    fClusterHistograms[i] = 0;
 }
 
 AliROCESDAnalysisSelector::~AliROCESDAnalysisSelector()
@@ -67,12 +64,6 @@ AliROCESDAnalysisSelector::~AliROCESDAnalysisSelector()
   //
   // Destructor
   //
-
-  if (fClusterHistograms) {
-    delete fClusterHistograms;
-    fClusterHistograms = 0;    
-  }
-
 }
 
 void AliROCESDAnalysisSelector::SlaveBegin(TTree* tree)
@@ -80,9 +71,6 @@ void AliROCESDAnalysisSelector::SlaveBegin(TTree* tree)
   //
   
   AliSelector::SlaveBegin(tree);
-
-  fClusterHistograms = new AliTPCClusterHistograms("test","test");
-
 } 
 
 void AliROCESDAnalysisSelector::Init(TTree *tree)
@@ -94,11 +82,20 @@ void AliROCESDAnalysisSelector::Init(TTree *tree)
   // Init() will be called many times when running with PROOF.
 
   AliSelector::Init(tree);
+  
+  printf("Init called %p\n", (void*) fESDfriend);
 
   // Set branch address
-  if (tree)
+  if (tree) 
+  {
     tree->SetBranchAddress("ESDfriend", &fESDfriend);
-    
+  
+    tree->SetBranchStatus("*", 0);
+    tree->SetBranchStatus("fTracks.*", 1);
+    tree->SetBranchStatus("fTimeStamp", 1);
+    //tree->SetBranchStatus("fTracks.fCalibContainer", 0);
+  }
+
   if (fESDfriend != 0)
     AliDebug(AliLog::kInfo, "INFO: Found ESDfriend branch in chain.");
 }
@@ -128,10 +125,64 @@ Bool_t AliROCESDAnalysisSelector::Process(Long64_t entry)
     return kFALSE;
   }
 
+  if (fESD->GetNumberOfTracks() != fESDfriend->GetNumberOfTracks())
+  {
+    AliDebug(AliLog::kError, Form("Event %lld: Number of tracks differ between ESD (%d) and ESDfriend (%d)! Skipping event!\n", entry, fESD->GetNumberOfTracks(), fESDfriend->GetNumberOfTracks()))
+    return kFALSE;
+  }  
+
   fESD->SetESDfriend(fESDfriend);
 
+  Int_t flag = ProcessEvent(entry, kFALSE);
+  if (flag == 1)
+    ProcessEvent(entry, kTRUE);
+
+  // TODO This should not be needed, the TTree::GetEntry() should take care of this, maybe because it has a reference member, to be analyzed
+  // if the ESDfriend is not deleted we get a major memory leak
+  // here the esdfriend seems to be also deleted, very weird behaviour....
+
+  delete fESD;
+  fESD = 0;
+  
+  //delete fESDfriend;
+  //fESDfriend = 0;
+
+  return kTRUE;
+}
+
+Int_t AliROCESDAnalysisSelector::ProcessEvent(Long64_t entry, Bool_t detailedHistogram)
+{
+  //
+  // Looping over tracks and clusters in event and filling histograms 
+  //
+  // - if detailedHistogram = kTRUE special histograms are saved (in fObjectsToSave)
+  // - the method returns 
+  //   1 : if a "flash" is detected something special in this event
+  //   
+
+  // save maximum 50 objects
+  if (detailedHistogram) 
+    if (fObjectsToSave->GetEntries() > 50) 
+      return 0;
+      
+  // for saving single events
+  AliTPCClusterHistograms* clusterHistograms[kTPCSectors];
+  for (Int_t i=0; i<kTPCSectors; i++) 
+    clusterHistograms[i] = 0;  
+
+  Bool_t intToReturn = 0;
+
   Int_t nTracks = fESD->GetNumberOfTracks();
   
+  Int_t nSkippedSeeds = 0;
+  Int_t nSkippedTracks = 0;
+
+  // for "flash" detection
+  Int_t nClusters = 0;
+  Float_t clusterQtotSumVsTime[250];  
+  for (Int_t z=0; z<250; z++)
+    clusterQtotSumVsTime[z] = 0;
+  
   // loop over esd tracks
   for (Int_t t=0; t<nTracks; t++)
   {
@@ -152,28 +203,148 @@ Bool_t AliROCESDAnalysisSelector::Process(Long64_t entry)
     const AliTPCseed* seed = dynamic_cast<const AliTPCseed*> (friendtrack->GetCalibObject(0));
     if (!seed)
     {
-      AliDebug(AliLog::kError, Form("ERROR: Could not retrieve seed of track %d.", t));
+      AliDebug(AliLog::kDebug, Form("ERROR: Could not retrieve seed of track %d.", t));
+      nSkippedSeeds++;
+      continue;
+    }
+    
+    if (!AcceptTrack(seed, fMinNumberOfRowsIsTrack))
+    {
+      AliDebug(AliLog::kDebug, Form("INFO: Rejected track %d.", t));
+      nSkippedTracks++;
       continue;
     }
     
     for (Int_t clusterID = 0; clusterID < 160; clusterID++)
     {
       AliTPCclusterMI* cluster = seed->GetClusterPointer(clusterID);
-      if (!cluster)
+      if (!cluster) 
+        continue;
+      
+      Int_t detector = cluster->GetDetector();
+      
+      if (detector < 0 || detector >= kTPCSectors) 
       {
-        //AliDebug(AliLog::kError, Form("ERROR: Could not retrieve cluster %d of track %d.", clusterID, t));
+        AliDebug(AliLog::kDebug, Form("We found a cluster from invalid sector %d", detector));
         continue;
       }
-      
-      //AliDebug(AliLog::kDebug, Form("We found a cluster from sector %d", cluster->GetDetector()));
 
-      if (cluster->GetDetector() != 5)
-        continue;
-      
-      fClusterHistograms->FillCluster(cluster);
+      if (!detailedHistogram) {
+        // TODO: find a clever way to handle the time      
+       Int_t time = 0;
+       
+       if (fESD->GetTimeStamp()>1160000000)
+         time = fESD->GetTimeStamp();      
+       
+       if (!fClusterHistograms[detector])
+         fClusterHistograms[detector] = new AliTPCClusterHistograms(detector,"",time,time+5*60*60);
+       
+       if (!fClusterHistograms[detector+kTPCSectors])
+         fClusterHistograms[detector+kTPCSectors] = new AliTPCClusterHistograms(detector,"",time,time+5*60*60, kTRUE);
+       
+       fClusterHistograms[detector]->FillCluster(cluster, time);
+       fClusterHistograms[detector+kTPCSectors]->FillCluster(cluster, time);
+       
+       Int_t z = Int_t(cluster->GetZ()); 
+       if (z>=0 && z<250) {
+         nClusters++;
+         clusterQtotSumVsTime[z] += cluster->GetQ();
+       }
+      } // end of if !detailedHistograms
+      else {
+       // if we need the detailed histograms for this event
+       if (!clusterHistograms[detector])
+         clusterHistograms[detector] = new AliTPCClusterHistograms(detector, Form("flash_entry%d", entry));
+       
+       clusterHistograms[detector]->FillCluster(cluster);
+      }
+    }
+    
+    for (Int_t i=0; i<kTPCHists; i++) 
+      if (fClusterHistograms[i]) 
+        fClusterHistograms[i]->FillTrack(seed);
+
+  }
+  
+  // check if there's a very large q deposit ("flash")
+  if (!detailedHistogram) {
+    for (Int_t z=0; z<250; z++) {
+      if (clusterQtotSumVsTime[z] > 150000) {
+       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()));
+               intToReturn = 1;
+      }
     }
   }
+  else {
+    for (Int_t i=0; i< kTPCSectors; i++) {
+      if (clusterHistograms[i]) {
+        fObjectsToSave->Add(clusterHistograms[i]);
+      }
+    }    
+  }
+
+//   if (nSkippedSeeds > 0)
+//     printf("WARNING: The seed was not found for %d out of %d tracks.\n", nSkippedSeeds, nTracks);
+//   if (nSkippedTracks > 0)
+//     printf("INFO: Rejected %d out of %d tracks.\n", nSkippedTracks, nTracks);
    
+  return intToReturn;
+}
+
+Bool_t AliROCESDAnalysisSelector::AcceptTrack(const AliTPCseed* track, Int_t minRowsIncluded) {
+  //
+  // check if the track should be accepted.
+  //
+  const Int_t   kMinClusters = 5;
+  const Float_t kMinRatio    = 0.75;
+  const Float_t kMax1pt      = 0.5;
+
+  Int_t  nRowsUsedByTracks = 0;
+  Bool_t rowIncluded[96];
+  
+  Float_t totalQtot = 0;
+  Int_t   nClusters = 0;
+
+  for(Int_t r=0; r<96; r++) 
+    rowIncluded[r] = kFALSE;
+  
+  for (Int_t clusterID = 0; clusterID < 160; clusterID++) {
+    AliTPCclusterMI* cluster = track->GetClusterPointer(clusterID);
+    
+    if (!cluster) 
+      continue;
+    
+    Float_t qTot =   cluster->GetQ();    
+    
+    nClusters++;
+    totalQtot += qTot;
+
+    if (!rowIncluded[cluster->GetRow()]) {
+      nRowsUsedByTracks++;
+      rowIncluded[cluster->GetRow()] = kTRUE;
+    }
+  }
+
+  Float_t meanQtot = totalQtot/nClusters;
+
+  if (meanQtot<70)
+    return kFALSE;
+
+  if (nRowsUsedByTracks < minRowsIncluded)
+    return kFALSE;
+  
+  //  printf(Form("    TRACK: n clusters = %d,  n pad rows = %d \n",track->GetNumberOfClusters(), nRowsUsedByTracks));
+  
+  if (track->GetNumberOfClusters()<kMinClusters) return kFALSE;
+  Float_t ratio = track->GetNumberOfClusters()/(track->GetNFoundable()+1.);
+  if (ratio<kMinRatio) return kFALSE;
+  Float_t mpt = track->Get1Pt();
+  if (TMath::Abs(mpt)>kMax1pt) return kFALSE;
+
+  //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
+  //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
+  //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
+  
   return kTRUE;
 }
 
@@ -181,15 +352,75 @@ void AliROCESDAnalysisSelector::SlaveTerminate()
 {
   //
   
-  fOutput->Add(fClusterHistograms);
+  if (fOutput)
+  {
+    for (Int_t i=0; i<kTPCHists; i++)
+      if (fClusterHistograms[i])
+        fOutput->Add(fClusterHistograms[i]);
+  }
 } 
 
 void AliROCESDAnalysisSelector::Terminate()
 {
+  // 
+  // read the objects from the output list and write them to a file
+  // the filename is modified by the object comment passed in the tree info or input list
   //
-    
-  TFile* file = TFile::Open("rocESD.root", "RECREATE");
+
+  if (fOutput)
+  {  
+    fOutput->Print();
+        
+    for (Int_t i=0; i<kTPCSectors; i++)
+      fClusterHistograms[i] = dynamic_cast<AliTPCClusterHistograms*> (fOutput->FindObject(AliTPCClusterHistograms::FormDetectorName(i, kFALSE)));
+    for (Int_t i=0; i<kTPCSectors; i++)
+      fClusterHistograms[kTPCSectors+i] = dynamic_cast<AliTPCClusterHistograms*> (fOutput->FindObject(AliTPCClusterHistograms::FormDetectorName(i, kTRUE)));
+  }
   
-  fClusterHistograms->SaveHistograms();
+  TNamed* comment = 0;
+  if (fTree && fTree->GetUserInfo())
+    comment = dynamic_cast<TNamed*>(fTree->GetUserInfo()->FindObject("comment"));
+  if (!comment && fInput)
+    comment = dynamic_cast<TNamed*>(fInput->FindObject("comment"));
+
+  if (comment)
+  {
+    AliDebug(AliLog::kInfo, Form("INFO: Found comment in input list: %s", comment->GetTitle()));
+  }
+  else
+    return;
+
+  TFile* file = TFile::Open(Form("rocESD_%s.root",comment->GetTitle()), "RECREATE");
+
+  for (Int_t i=0; i<kTPCHists; i++)
+    if (fClusterHistograms[i]) {
+      fClusterHistograms[i]->SaveHistograms();
+      TCanvas* c = fClusterHistograms[i]->DrawHistograms();
+                       TString dir;
+                       dir.Form("WWW/%s/%s", comment->GetTitle(), c->GetName());
+                       gSystem->mkdir(dir, kTRUE);
+      c->SaveAs(Form("%s/plots_%s_%s.eps",dir.Data(),comment->GetTitle(),c->GetName()));
+      c->SaveAs(Form("%s/plots_%s_%s.gif",dir.Data(),comment->GetTitle(),c->GetName()));
+
+      c->Close();
+      delete c;
+    }
+
+  gDirectory->mkdir("saved_objects");
+  gDirectory->cd("saved_objects");
+
+  for (Int_t i=0; i<fObjectsToSave->GetSize(); i++) {
+    if (fObjectsToSave->At(i)) {
+      AliTPCClusterHistograms* clusterHistograms = dynamic_cast<AliTPCClusterHistograms*> (fObjectsToSave->At(i));
+      if (clusterHistograms)
+       clusterHistograms->SaveHistograms();
+      else
+       fObjectsToSave->At(i)->Write();
+    }
+  }
+
+  gDirectory->cd("../");
+
+
   file->Close();
-} 
+}