HMPID Analysis task improved
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDAnalysisTask.cxx
index 8f1d768..888f4b3 100644 (file)
 #ifndef AliHMPIDAnalysisTASK_CXX
 #define AliHMPIDAnalysisTASK_CXX
 
-#include "AliHMPIDAnalysisTask.h"
-#include "TCanvas.h"
-#include "AliStack.h"
-#include "TParticle.h"
-#include "Riostream.h"
-#include "TH1I.h"
-#include "TH2F.h"
-#include "TH1F.h"
-#include "AliMCEvent.h"
+
+#include "TH1.h"
+#include "TH2.h"
 #include "AliAnalysisManager.h"
-#include "AliESDEvent.h"
-#include "TChain.h"
+#include "AliESDInputHandler.h"
 #include "AliESDtrack.h"
 #include "AliLog.h"
+#include "AliHMPIDAnalysisTask.h"
 
 ClassImp(AliHMPIDAnalysisTask)
 
 //__________________________________________________________________________
 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask() :
-  fHistList(0x0),
-  fHistEventsProcessed(0x0),
+  fESD(0x0),fHmpHistList(0x0),
   fNevts(0),
   fTrigNevts(0),
-  fTrigger(0)
+  fTrigger(0),
+  fHmpInner(0x0),fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
+  fHmpMipTrkDist(0x0),fHmpMipTrkDistX(0x0),fHmpMipTrkDistY(0x0),fHmpMipCharge3cm(0x0),fHmpMipCharge1cm(0x0),fHmpNumPhots(0x0),
+  fHmpTrkFlags(0x0)
 {
   //
   //Default ctor
   //
+  
 }
 //___________________________________________________________________________
 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
   AliAnalysisTaskSE(name),
-  fHistList(0x0),
-  fHistEventsProcessed(0x0),
-  fNevts(0),
+  fESD(0),fHmpHistList(0x0), fNevts(0),
   fTrigNevts(0),
-  fTrigger(0)
+  fTrigger(0),
+  fHmpInner(0x0),fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
+  fHmpMipTrkDist(0x0),fHmpMipTrkDistX(0x0),fHmpMipTrkDistY(0x0),fHmpMipCharge3cm(0x0),fHmpMipCharge1cm(0x0),fHmpNumPhots(0x0),
+  fHmpTrkFlags(0x0)
 {
   //
   // Constructor. Initialization of Inputs and Outputs
   //
-  Info("AliHMPIDAnalysisTask","Calling Constructor");
-
-  /*
-    DefineInput(0) and DefineOutput(0)
-    are taken care of by AliAnalysisTaskSE constructor
-  */
-  DefineOutput(1,TH1I::Class());
-  DefineOutput(2,TList::Class());
-}
-
-//___________________________________________________________________________
-AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c) 
-{
-  //
-  // Assignment operator
-  //
-  if (this!=&c) {
-    AliAnalysisTaskSE::operator=(c) ;
-    fHistList = c.fHistList ;
-    fHistEventsProcessed = c.fHistEventsProcessed;
-    fNevts = c.fNevts;
-    fTrigNevts = c.fTrigNevts;
-    fTrigger = c.fTrigger;
-  }
-  return *this;
-}
-
-//___________________________________________________________________________
-AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
-  AliAnalysisTaskSE(c),
-  fHistList(c.fHistList),
-  fHistEventsProcessed(c.fHistEventsProcessed),
-  fNevts(c.fNevts),
-  fTrigNevts(c.fTrigNevts),
-  fTrigger(c.fTrigger)
-{
-  //
-  // Copy Constructor
-  //
-}
-
+  
+  DefineOutput(0,TList::Class());
+ }
 //___________________________________________________________________________
 AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
   //
   //destructor
   //
   Info("~AliHMPIDAnalysisask","Calling Destructor");
-  if (fHistEventsProcessed) delete fHistEventsProcessed ;
-  if (fHistList) {fHistList->Clear(); delete fHistList;}
+  if (fHmpHistList) {fHmpHistList->Clear(); delete fHmpHistList;}
 }
-
+//___________________________________________________________________________
+void AliHMPIDAnalysisTask::ConnectInputData(Option_t *) 
+{
+  // Connect ESD here
+  
+  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+  if (!esdH) {
+      AliDebug(2,Form("ERROR: Could not get ESDInputHandler")); 
+    } else
+      fESD = esdH->GetEvent();
+ }
 //_________________________________________________
-void AliHMPIDAnalysisTask::UserExec(Option_t *)
+void AliHMPIDAnalysisTask::Exec(Option_t *)
 {
   //
   // Main loop function, executed on Event basis
   //
-  Info("UserExec","") ;
-
-  AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
-  if (!fESD) {
-    Error("UserExec","NO ESD FOUND!");
-    return;
-  }
-  fNevts++;
-  if(fESD->GetTriggerMask()&fTrigger == fTrigger) fTrigNevts++;
-
-
   for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
 
     AliESDtrack* track = fESD->GetTrack(iTrack);
@@ -133,47 +90,52 @@ void AliHMPIDAnalysisTask::UserExec(Option_t *)
     Double_t rin[3], rout[3];
     track->GetInnerXYZ(rin);
     track->GetOuterXYZ(rout);
+        /*
+     ((TH2F *)fHmpHistList->At(22))->Fill(rin[0],rin[1]);
+     ((TH2F *)fHmpHistList->At(23))->Fill(rout[0],rout[1]);
+*/
+
+     if(track->GetHMPIDsignal() == -20) fHmpTrkFlags->Fill(0);
+     else if(track->GetHMPIDsignal() == -9) fHmpTrkFlags->Fill(1);
+     else if(track->GetHMPIDsignal() == -5) fHmpTrkFlags->Fill(2);
+     else if(track->GetHMPIDsignal() == -11) fHmpTrkFlags->Fill(3);
+     else if(track->GetHMPIDsignal() == -22) fHmpTrkFlags->Fill(4);
+     else fHmpTrkFlags->Fill(4);
   
-     ((TH2F *)fHistList->At(22))->Fill(rin[0],rin[1]);
-     ((TH2F *)fHistList->At(23))->Fill(rout[0],rout[1]);
-
-
-     TH1F *h = ((TH1F*)(fHistList->At(0)));
-     if(track->GetHMPIDsignal() == -20) h->Fill(0);
-     else if(track->GetHMPIDsignal() == -9) h->Fill(1);
-     else if(track->GetHMPIDsignal() == -5) h->Fill(2);
-     else if(track->GetHMPIDsignal() == -11) h->Fill(3);
-     else if(track->GetHMPIDsignal() == -22) h->Fill(4);
-     else h->Fill(4);
-    
-     if(track->GetHMPIDcluIdx() == 99099999) continue;
-     if(track->GetHMPIDcluIdx() < 0) continue;
-
+    if(track->GetHMPIDsignal()== -20) continue;
+    if(track->GetHMPIDcluIdx() < 0) continue;
+   
     Int_t ch = track->GetHMPIDcluIdx()/1000000; 
-
     Float_t x, y; Int_t q, nph; 
-    track->GetHMPIDmip(x,y,q,nph);
-
-    if(x==0 && y==0) continue;
     Float_t xpc, ypc, th, ph;
-    track->GetHMPIDtrk(xpc,ypc,th,ph); //special setting in local cosmic rec!!!
-                                       //do not use with standard rec
-
+    
+    track->GetHMPIDmip(x,y,q,nph);
+    track->GetHMPIDtrk(xpc,ypc,th,ph);
+     
+    if(x==0 && y==0 && xpc == 0 && ypc == 0) continue;
+    //Printf("%s  %s Good track is found",(char*)__FILE__,__LINE__);
+     
     Double_t dist = TMath::Sqrt( (xpc-x)*(xpc-x) + (ypc - y)*(ypc - y));
-
-     if(ch >=0 && ch < 7) {
-     if(dist < 3) ((TH2F *)fHistList->At(ch+1))->Fill(x,y);
-     ((TH1F *)fHistList->At(ch+7+1))->Fill(q);
-     ((TH2F *)fHistList->At(ch+14+1))->Fill(dist,q);
-     }
-  
+    fHmpMipTrkDist->Fill(dist);
+    fHmpMipTrkDistX->Fill(xpc-x);
+    fHmpMipTrkDistY->Fill(ypc - y);
+    if(dist<=3.0) fHmpMipCharge3cm->Fill(q);
+    
+    if(track->GetHMPIDsignal() > 0 )
+    {
+      Printf("EvtNumInFile: %d HMPID ThetaC: %lf x: %lf xpc: %lf y: %lf ypx: %lf Q: %d nPh: %d IncTheta; %lf IncPhi: %lf",fESD->GetEventNumberInFile(),track->GetHMPIDsignal(),x,xpc,y,ypc,q,nph,th,ph);
+      Double_t pHmp[3] = {0},pmod = 0;if(track->GetOuterHmpPxPyPz(pHmp))  pmod = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]); 
+      fHmpPesdPhmp->Fill(track->P(),pmod);  
+      if(dist<=1.0)  fHmpMipCharge1cm->Fill(q);
+      fHmpNumPhots->Fill(nph);
+      fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
+      fHmpCkovPesd->Fill(pmod,track->GetHMPIDsignal());
+    }//there is signal
+    
     }//track loop
-
-  fHistEventsProcessed->Fill(0);
   
   /* PostData(0) is taken care of by AliAnalysisTaskSE */
-  PostData(1,fHistEventsProcessed) ;
-  PostData(2,fHistList) ;
+  PostData(0,fHmpHistList) ;
 }
 
 
@@ -191,25 +153,53 @@ void AliHMPIDAnalysisTask::Terminate(Option_t*)
 
 
 //___________________________________________________________________________
-void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
+void AliHMPIDAnalysisTask::CreateOutputObjects() {
   //
   //HERE ONE CAN CREATE OUTPUT OBJECTS
   //TO BE SET BEFORE THE EXECUTION OF THE TASK
   //
-  Info("CreateOutputObjects","CreateOutputObjects of task %s", GetName());
 
   //slot #1
-  OpenFile(1);
-  fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
-
-   OpenFile(2);
-   fHistList = new TList();
+   OpenFile(0);
+   fHmpHistList = new TList();
+   fHmpInner =new TH2F("fHmpInner","HMPID: Inner track XY;X (cm);Y(cm)",800,-400,400,800,-400,400);
+   fHmpHistList->Add(fHmpInner);
+   
+   fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
+   fHmpHistList->Add(fHmpPesdPhmp);
+   
+   fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
+   fHmpHistList->Add(fHmpCkovPesd);
+   
+   fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
+   fHmpHistList->Add(fHmpCkovPhmp);  
    
+   fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",800,-20,20);
+   fHmpHistList->Add(fHmpMipTrkDist);
+   fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
+   fHmpHistList->Add(fHmpMipTrkDistX);
+   fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
+   fHmpHistList->Add(fHmpMipTrkDistY);
+   
+   fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
+   fHmpHistList->Add(fHmpMipCharge3cm);
+   
+   fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
+   fHmpHistList->Add(fHmpMipCharge1cm);
+   
+   fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
+   fHmpHistList->Add(fHmpNumPhots);
+   
+   fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
+   TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
+   for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
+   fHmpHistList->Add(fHmpTrkFlags);
+   /*
    //0
    TH1F *trkH = new TH1F("trkH","signal flags in HMPID",6,0,6);
    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
    for(Int_t ibin = 0; ibin < 6; ibin++) trkH->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
-   fHistList->Add(trkH);
+   fHmpHistList->Add(trkH);
 
   
    TH2F *mod[7], *dq[7];
@@ -220,33 +210,33 @@ void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
    mod[i] = new TH2F(Form("mod%i",i),Form("MIP position in chamber %i",i),180,0,180,180,0,180);
    mod[i]->SetMarkerStyle(8);
    mod[i]->SetMarkerColor(2);
-   fHistList->Add(mod[i]);
+   fHmpHistList->Add(mod[i]);
    }
    //8-14
    for(Int_t i=0; i< 7 ; i++) {//to reserve the histo sorting
    q[i] = new TH1F(Form("q%i",i),Form("MIP charge distribution in chamber %i",i),5000,0,5000);
    q[i]->SetMarkerColor(2);
-   fHistList->Add(q[i]);
+   fHmpHistList->Add(q[i]);
    }
    //15-21
    for(Int_t i=0; i< 7 ; i++) {//to reserve the histo sorting
    dq[i] = new TH2F(Form("dq%i",i),Form("#Delta(mip-track) vs Q_{mip} in chamber %i",i),1000,0,100,5000,0,5000),
    dq[i]->SetMarkerStyle(6);
-   fHistList->Add(dq[i]);
+   fHmpHistList->Add(dq[i]);
    }
    //22
    TH2F *inner = new TH2F("inner","inner track XY",800,-400,400,800,-400,400);
    inner->SetMarkerStyle(6);
    inner->SetXTitle("X cm");
    inner->SetYTitle("Y cm");
-   fHistList->Add(inner); 
+   fHmpHistList->Add(inner); 
    //23
    TH2F *outer = new TH2F("outer","outer track XY",800,-400,400,800,-400,400);
    outer->SetMarkerStyle(6);
    outer->SetXTitle("X cm");
    outer->SetYTitle("Y cm");
-   fHistList->Add(outer);
-
+   fHmpHistList->Add(outer);
+*/
 }
 
 #endif