PWG0/multPbPb/AliAnalysisTaskTriggerStudy.cxx

   1 // AliAnalysisTaskTriggerStudy
   2
   3 // Author: Michele Floris, CERN
   4 // TODO:
   5 // - Add chi2/cluster plot for primary, secondaries and fakes
   6
   7
   8 #include "AliAnalysisTaskTriggerStudy.h"
   9 #include "AliESDInputHandler.h"
  10 #include "AliHistoListWrapper.h"
  11 #include "AliAnalysisManager.h"
  12 #include "AliMCEvent.h"
  13 #include "AliStack.h"
  14 #include "TH1I.h"
  15 #include "TH3D.h"
  16 #include "AliMCParticle.h"
  17 #include "AliGenEventHeader.h"
  18 #include "AliESDCentrality.h"
  19
  20 #include <iostream>
  21 #include "AliTriggerAnalysis.h"
  22 #include "AliMultiplicity.h"
  23 #include "TFile.h"
  24 #include "AliLog.h"
  25 #include "AliESDtrackCuts.h"
  26 #include "AliESDVZERO.h"
  27 #include "TH2F.h"
  28
  29 using namespace std;
  30
  31 ClassImp(AliAnalysisTaskTriggerStudy)
  32
  33 //const char * AliAnalysisTaskTriggerStudy::kVDNames[] = {"C0SM1","C0SM2","C0VBA","C0VBC","C0OM2"};
  34 const char * AliAnalysisTaskTriggerStudy::kVDNames[] = {"C0SM1","C0SM2","C0VBA","C0VBC"};//,"C0OM2"};
  35
  36 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy()
  37 : AliAnalysisTaskSE("TaskTriggerStudy"),
  38   fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
  39 {
  40   // constructor
  41
  42   DefineOutput(1, AliHistoListWrapper::Class());
  43
  44 }
  45 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const char * name)
  46   : AliAnalysisTaskSE(name),
  47     fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
  48 {
  49   //
  50   // Standard constructur which should be used
  51   //
  52
  53   DefineOutput(1, AliHistoListWrapper::Class());
  54
  55 }
  56
  57 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const AliAnalysisTaskTriggerStudy& obj) :
  58   AliAnalysisTaskSE(obj) ,fESD (0), fIsMC(0), fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
  59 {
  60   //copy ctor
  61   fESD = obj.fESD ;
  62   fHistoList = obj.fHistoList;
  63   fTriggerAnalysis = obj.fTriggerAnalysis;
  64   fHistoSuffix = obj.fHistoSuffix;
  65   fNTrackletsCut = obj.fNTrackletsCut;
  66   fNTrackletsCutKine = obj.fNTrackletsCutKine;
  67   fRejectBGWithV0 = obj.fRejectBGWithV0;
  68 }
  69
  70 AliAnalysisTaskTriggerStudy::~AliAnalysisTaskTriggerStudy(){
  71   // destructor
  72
  73   if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
  74     if(fHistoList) {
  75       delete fHistoList;
  76       fHistoList = 0;
  77     }
  78     if(fTriggerAnalysis) {
  79       delete fTriggerAnalysis;
  80       fHistoList = 0;
  81     }
  82   }
  83   // Histo list should not be destroyed: fListWrapper is owner!
  84
  85 }
  86 void AliAnalysisTaskTriggerStudy::UserCreateOutputObjects()
  87 {
  88   // Called once
  89   fHistoList = new AliHistoListWrapper("histoList","histogram list for trigger studies");
  90   fTriggerAnalysis = new AliTriggerAnalysis();
  91 }
  92
  93
  94 void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
  95 {
  96   // User code
  97
  98   // FIXME: make sure you have the right cuts here
  99
 100   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 101   PostData(1,fHistoList);
 102
 103   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 104   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
 105     AliFatal("Not processing ESDs");
 106   }
 107
 108   // get the multiplicity object
 109   const AliMultiplicity* mult = fESD->GetMultiplicity();
 110   Int_t ntracklets = mult->GetNumberOfTracklets();
 111   // Get Number of tracks
 112   Int_t ntracks    = AliESDtrackCuts::GetReferenceMultiplicity(fESD,kTRUE); // tpc only
 113
 114   // Get V0 Multiplicity
 115   AliESDVZERO* esdV0 = fESD->GetVZEROData();
 116   Float_t multV0A=esdV0->GetMTotV0A();
 117   Float_t multV0C=esdV0->GetMTotV0C();
 118   Float_t multV0 = multV0A+multV0C;
 119
 120   // Get number of clusters in layer 1
 121   Float_t outerLayerSPD = mult->GetNumberOfITSClusters(1);
 122   Float_t innerLayerSPD = mult->GetNumberOfITSClusters(0);
 123   Float_t totalClusSPD = outerLayerSPD+innerLayerSPD;
 124
 125   GetHistoSPD1  ("All", "All events before any selection")->Fill(outerLayerSPD);
 126
 127
 128   // Physics selection
 129   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 130   if(!isSelected) return;
 131   // FIXME trigger classes
 132   // if ( !(fESD->IsTriggerClassFired("CMBS2A-B-NOPF-ALL")|| fESD->IsTriggerClassFired("CMBS2C-B-NOPF-ALL") || fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL")) ) return;
 133
 134   GetHistoSPD1  ("AllPSNoPrino", "All events after physsel, before Francesco")->Fill(outerLayerSPD);
 135
 136
 137   // Francesco's cuts
 138   const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC();
 139   if(vtxESDTPC->GetNContributors()<1) return;
 140   if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0)))     return;
 141   const AliESDVertex * vtxESDSPD= fESD->GetPrimaryVertexSPD();
 142   Float_t tpcContr=vtxESDTPC->GetNContributors();
 143
 144
 145
 146   // GetT0 Stuff
 147   const Double32_t *meanT0   = fESD->GetT0TOF();
 148   const Double32_t  meanT0A  = 0.001* meanT0[1];
 149   const Double32_t  meanT0C  = 0.001* meanT0[2];
 150   const Double32_t  meanT0AC = 0.001* meanT0[0];
 151   Double32_t  T0Vertex = fESD->GetT0zVertex();
 152   //  Double32_t  *ampT0 =Esdevent ->GetT0amplitude();
 153
 154 //   cut1yesF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && francescoscut
 155 // cut1notF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && ! francescoscut
 156 // cut2 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) <=-0.7) && meanC > -2) )
 157 // cut3 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) < 0.7 && (meanC-meanA) > -0.7 ) )
 158 //   cut4 = ( (meanC < 95. && meanA < 95.) && (meanA < -2.)
 159
 160   Bool_t cut1T0 =  ( (meanT0C < 95. && meanT0A < 95.) && (meanT0C < -2.) );
 161   Bool_t cut2T0 = ( (meanT0C < 95. && meanT0A < 95.) && ( (meanT0C-meanT0A) <=-0.7) && meanT0C > -2) ;
 162
 163   if(ntracklets > fNTrackletsCut) return;
 164
 165   // Reset histo suffix and fill reference histograms without any suffix
 166   fHistoSuffix = "";
 167
 168   // Fast or in the outer layer
 169   Int_t nFastOrOnline  = fTriggerAnalysis->SPDFiredChips(fESD, 1, 0, 2); // offline
 170   Int_t nFastOrOffline = fTriggerAnalysis->SPDFiredChips(fESD, 0, 0, 2); // online
 171
 172   Bool_t c0sm1 = nFastOrOffline >= 1;
 173   Bool_t c0sm2 = nFastOrOffline >= 2;
 174   Bool_t c0sm3 = nFastOrOffline >= 3;
 175   Bool_t c0sm4 = nFastOrOffline >= 4;
 176   Bool_t c0sm5 = nFastOrOffline >= 5;
 177
 178   // V0 triggers
 179   Bool_t c0v0A       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0A);
 180   Bool_t c0v0C       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0C);
 181   Bool_t v0AHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 182   Bool_t v0CHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 183
 184   Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0ABG);
 185   Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0CBG);
 186   Bool_t v0BG  = v0ABG || v0CBG;
 187
 188
 189   // TOF triggers
 190   // FIXME: move to triggeranalysis?
 191   AliESDHeader*h = fESD->GetHeader(); // taken the header from AliESDEvent
 192   Bool_t c0OM2 = h->IsTriggerInputFired("0OM2"); // thr >= 2 (input 19)
 193   Bool_t c0OM3 = h->IsTriggerInputFired("0OM3"); // thr >= 3 (input 20)
 194
 195   // ZDC triggers
 196   Bool_t zdcA   = kFALSE;
 197   Bool_t zdcC   = kFALSE;
 198   Bool_t zdcBar = kFALSE;
 199
 200   if (!fIsMC) {
 201     // If it's data, we use the TDCs
 202     zdcA   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kASide, kTRUE, kFALSE) ;
 203     zdcC   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCSide, kTRUE, kFALSE) ;
 204     zdcBar = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCentralBarrel) ;
 205   } else {
 206     // If it's MC, we use the energy
 207     Double_t minEnergy = 0;
 208     AliESDZDC *esdZDC = fESD->GetESDZDC();
 209     Double_t  zNCEnergy = esdZDC->GetZDCN1Energy();
 210     Double_t  zPCEnergy = esdZDC->GetZDCP1Energy();
 211     Double_t  zNAEnergy = esdZDC->GetZDCN2Energy();
 212     Double_t  zPAEnergy = esdZDC->GetZDCP2Energy();
 213     // zdcA = (zNAEnergy>minEnergy || zPAEnergy>minEnergy);
 214     // zdcC = (zNCEnergy>minEnergy || zPCEnergy>minEnergy);
 215     zdcA = (zNAEnergy>minEnergy);
 216     zdcC = (zNCEnergy>minEnergy);
 217   }
 218
 219
 220   // Some macros for the online triggers
 221   Bool_t cMBS2A = c0sm2 && c0v0A;
 222   Bool_t cMBS2C = c0sm2 && c0v0C;
 223   Bool_t cMBAC  = c0v0A && c0v0C;
 224
 225
 226   Bool_t vdArray[kNVDEntries];
 227   vdArray[kVDC0MBS1] = c0sm1;
 228   vdArray[kVDC0MBS2] = c0sm2;
 229   vdArray[kVDC0VBA]  = c0v0A;
 230   vdArray[kVDC0VBC]  = c0v0C;
 231   //vdArray[kVDC0OM2]  = c0OM2;
 232
 233   // Plots requested by Jurgen on 18/11/2010
 234   // FIXME: will skip everything else
 235
 236   GetHistoSPD1  ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(outerLayerSPD);
 237   GetHistoTracks("PhysSel", "All events after physics selection and Francesco's cut")->Fill(ntracks);
 238   GetHistoV0M   ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(multV0);
 239   if(c0v0A && c0v0C) {
 240     GetHistoSPD1  ("V0AND", "V0A & V0C")->Fill(outerLayerSPD);
 241     GetHistoTracks("V0AND", "V0A & V0C")->Fill(ntracks);
 242     GetHistoV0M   ("V0AND", "V0A & V0C")->Fill(multV0);
 243   }
 244   if((c0v0A && !c0v0C) || (!c0v0A && c0v0C)) {
 245     GetHistoSPD1  ("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(outerLayerSPD);
 246     GetHistoTracks("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(ntracks);
 247     GetHistoV0M   ("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(multV0);
 248   }
 249   if(zdcA && zdcC) {
 250     GetHistoSPD1  ("ZDCAND", "ZDCA & ZDCC")->Fill(outerLayerSPD);
 251     GetHistoTracks("ZDCAND", "ZDCA & ZDCC")->Fill(ntracks);
 252     GetHistoV0M   ("ZDCAND", "ZDCA & ZDCC")->Fill(multV0);
 253   }
 254   if((c0v0A && c0v0C) && !(zdcA && zdcC)) {
 255     GetHistoSPD1  ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(outerLayerSPD);
 256     GetHistoTracks("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(ntracks);
 257     GetHistoV0M   ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(multV0);
 258   }
 259   if((c0v0A && c0v0C) && !zdcA && !zdcC) {
 260     GetHistoSPD1  ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(outerLayerSPD);
 261     GetHistoTracks("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(ntracks);
 262     GetHistoV0M   ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(multV0);
 263 }
 264   if((c0v0A && c0v0C) && ((!zdcA && zdcC) || (zdcA && !zdcC))) {
 265     GetHistoSPD1  ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(outerLayerSPD);
 266     GetHistoTracks("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(ntracks);
 267     GetHistoV0M   ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(multV0);
 268   }
 269
 270   if(((c0v0A && !c0v0C) || (!c0v0A && c0v0C)) && (zdcA && zdcC)) {
 271     GetHistoSPD1  ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(outerLayerSPD);
 272     GetHistoTracks("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(ntracks);
 273     GetHistoV0M   ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(multV0);
 274   }
 275
 276   if((c0v0A && c0v0C) && (zdcA && zdcC)) {
 277     GetHistoSPD1  ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(outerLayerSPD);
 278     GetHistoTracks("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(ntracks);
 279     GetHistoV0M   ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(multV0);
 280   }
 281
 282   if((c0v0A && zdcA && !zdcC && !c0v0C) || (c0v0C && zdcC && !zdcA && !c0v0A)) {
 283     GetHistoSPD1  ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(outerLayerSPD);
 284     GetHistoTracks("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(ntracks);
 285     GetHistoV0M   ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(multV0);
 286   }
 287
 288   // GetHistoCorrelationSPDTPCVz("All", "After physics selection and Francesco's cut")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
 289   // if(cut1T0)   GetHistoCorrelationSPDTPCVz("Cut1T0", "T0 Cut 1")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
 290   // if(cut2T0)   GetHistoCorrelationSPDTPCVz("Cut2T0", "T0 Cut 2")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
 291
 292   // GetHistoCorrelationContrTPCSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,tpcContr);
 293   // if(cut1T0)   GetHistoCorrelationContrTPCSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,tpcContr);
 294   // if(cut2T0)   GetHistoCorrelationContrTPCSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,tpcContr);
 295
 296   // GetHistoCorrelationTrackletsSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,ntracklets);
 297   // if(cut1T0)   GetHistoCorrelationTrackletsSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,ntracklets);
 298   // if(cut2T0)   GetHistoCorrelationTrackletsSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,ntracklets);
 299
 300
 301   return; // FIXME
 302
 303
 304   // Reject background
 305   if (v0BG && fRejectBGWithV0) {
 306     cout << "Rejection BG" << endl;
 307
 308     return;
 309   }
 310   // Fill global histos
 311   GetHistoTracklets("all","All events")->Fill(ntracklets);
 312   FillTriggerOverlaps("All", "All Events",vdArray);
 313
 314   // Fill some combination of trigger classes
 315   Bool_t cmbs1aOnline = fESD->IsTriggerClassFired("CMBS1A-B-NOPF-ALL");
 316   Bool_t cmbs1cOnline = fESD->IsTriggerClassFired("CMBS1C-B-NOPF-ALL");
 317   Bool_t cmbacOnline  = fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL");
 318
 319   Bool_t twoOutOfThree = kFALSE;
 320   if (cmbs1aOnline || cmbs1cOnline ||cmbacOnline) twoOutOfThree = kTRUE;
 321   if (twoOutOfThree) GetHistoTracklets("All_TwoOutOfThree" ,"Events 2-out-of-3 online" )->Fill(ntracklets);
 322
 323
 324   // loop over trigger classes in the event
 325   TObjArray * tokens = 0;
 326   if(fIsMC) {
 327     // in case of montecarlo I override the trigger class
 328     tokens = new TObjArray;
 329     tokens->SetOwner();
 330     //    tokens->Add(new TObjString("CINT1B-ABCE-NOPF-ALL"));
 331     tokens->Add(new TObjString("MC"));
 332   }
 333   else {
 334     TString trgClasses = fESD->GetFiredTriggerClasses();
 335     tokens = trgClasses.Tokenize("  ");
 336   }
 337   TIterator  * iter = (TIterator*) tokens->MakeIterator();
 338
 339   TString classes = fESD->GetFiredTriggerClasses();
 340
 341   // if (classes.Contains("SMH")) {
 342   //    tokens->Print();
 343   //    cout << classes.Data() << endl;
 344   // }
 345   //  iter->Reset();
 346   //Int_t itoken = 0;
 347   TObjString * tok=0;
 348   while((tok = (TObjString*) iter->Next())){
 349     // clean up trigger name
 350     TString trg = tok->GetString();
 351     trg.Strip(TString::kTrailing, ' ');
 352     trg.Strip(TString::kLeading, ' ');
 353
 354     fHistoSuffix = "_";
 355     fHistoSuffix += trg;
 356
 357     //    cout << itoken++ << " " << trg.Data() << endl;
 358     // continue;
 359     //    if (trg.Contains("SMH")) cout << itoken++ << " " << trg.Data() << endl;
 360
 361     // Fill tracklets
 362     GetHistoTracklets("All" ,"Events no offline trigger" )->Fill(ntracklets);
 363
 364
 365     // Fill histograms mismatchs
 366     // TODO: check mismatch trigger class
 367     if(nFastOrOffline != nFastOrOnline) {
 368       GetHistoTracklets("mismatchingFastOr", "Events where fast or offline differs from fast-or online")->Fill(ntracklets);
 369     }
 370
 371     if (c0v0A != v0AHW){
 372       GetHistoTracklets("mismatchingV0A", "Events where V0A offline differs from V0A online")->Fill(ntracklets);
 373     }
 374
 375     if (c0v0C != v0CHW){
 376       GetHistoTracklets("mismatchingV0C", "Events where V0C offline differs from V0C online")->Fill(ntracklets);
 377     }
 378
 379     // Fill a tracklet histo for each trigger type
 380     if(c0sm1)  GetHistoTracklets("c0sm1" ,"Events were trigger c0sm1 fired" )->Fill(ntracklets);
 381     if(c0sm2)  GetHistoTracklets("c0sm2" ,"Events were trigger c0sm2 fired" )->Fill(ntracklets);
 382     if(c0sm3)  GetHistoTracklets("c0sm3" ,"Events were trigger c0sm3 fired" )->Fill(ntracklets);
 383     if(c0sm4)  GetHistoTracklets("c0sm4" ,"Events were trigger c0sm4 fired" )->Fill(ntracklets);
 384     if(c0sm5)  GetHistoTracklets("c0sm5" ,"Events were trigger c0sm5 fired" )->Fill(ntracklets);
 385     if(c0OM2)  GetHistoTracklets("c0OM2" ,"Events were trigger c0OM2 fired" )->Fill(ntracklets);
 386     if(c0OM3)  GetHistoTracklets("c0OM3" ,"Events were trigger c0OM3 fired" )->Fill(ntracklets);
 387     if(c0v0A)  GetHistoTracklets("c0v0A" ,"Events were trigger c0v0A fired" )->Fill(ntracklets);
 388     if(c0v0C)  GetHistoTracklets("c0v0C" ,"Events were trigger c0v0C fired" )->Fill(ntracklets);
 389     if(cMBS2A) GetHistoTracklets("cMBS2A","Events were trigger cMBS2A fired")->Fill(ntracklets);
 390     if(cMBS2C) GetHistoTracklets("cMBS2C","Events were trigger cMBS2C fired")->Fill(ntracklets);
 391     if(cMBAC ) GetHistoTracklets("cMBAC", "Events were trigger cMBAC  fired")->Fill(ntracklets);
 392     if(zdcA )  GetHistoTracklets("cZDCA", "Events were trigger cZDCA  fired")->Fill(ntracklets);
 393     if(zdcC )  GetHistoTracklets("cZDCC", "Events were trigger cZDCC  fired")->Fill(ntracklets);
 394     if(!zdcA && !zdcC )  GetHistoTracklets("NoZDC", "Events were zdc trigger don't  fired")->Fill(ntracklets);
 395     if( (zdcA && !zdcC) || (!zdcA && zdcC) )  GetHistoTracklets("OneSideZDC", "Events were only 1 ZDC  fired")->Fill(ntracklets);
 396     if( (zdcA && zdcC) )  GetHistoTracklets("TwoSideZDC", "Events were both  ZDC  fired")->Fill(ntracklets);
 397
 398     if(twoOutOfThree) {
 399       if(!zdcA && !zdcC )  GetHistoTracklets("NoZDC_TwoOutOfThree", "Events were zdc trigger don't  fired")->Fill(ntracklets);
 400       if( (zdcA && !zdcC) || (!zdcA && zdcC) )  GetHistoTracklets("OneSideZDC_TwoOutOfThree", "Events were only 1 ZDC  fired")->Fill(ntracklets);
 401       if( (zdcA && zdcC) )  GetHistoTracklets("TwoSideZDC_TwoOutOfThree", "Events were both  ZDC  fired")->Fill(ntracklets);
 402     }
 403     if(zdcBar) GetHistoTracklets("cZDCBar","Events were trigger cZDCB  fired")->Fill(ntracklets);
 404     //  if() GetHistoTracklets("","Events were trigger  fired");
 405
 406     // Fill trigger overlaps
 407     FillTriggerOverlaps("All", "All Events in trigger class",vdArray);
 408
 409
 410     // Fill kinematic variables for peripheral events
 411     static AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(); // FIXME: change this ?
 412     if (ntracklets < fNTrackletsCutKine) {
 413       // 1. Loop on tracklets
 414       for(Int_t itracklet = 0; itracklet < ntracklets; itracklet++){
 415         GetHistoEta("All", Form("Tracklet #eta distribution, ntracklets < %d",fNTrackletsCutKine))->Fill(mult->GetEta(itracklet));
 416       }
 417       // 2. loop on tracks
 418       TObjArray * goodTracks = cuts->GetAcceptedTracks(fESD);
 419       TIterator * iterTracks = goodTracks->MakeIterator();
 420       AliESDtrack * esdTrack = 0;
 421       while ((esdTrack = (AliESDtrack*) iterTracks->Next())) {
 422         GetHistoPt("All", Form("Tracklet p_{T} distribution, ntracklets < %d",fNTrackletsCutKine))->Fill(esdTrack->Pt());
 423       }
 424       delete goodTracks;
 425     }
 426
 427
 428   }
 429   delete tokens;
 430
 431     // if (fIsMC) {
 432
 433
 434   //   if (!fMCEvent) {
 435   //     AliError("No MC info found");
 436   //   } else {
 437
 438   //     //loop on the MC event
 439   //     //      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
 440   //     Int_t offset    = fMCEvent->GetPrimaryOffset();
 441   //     Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries()+offset;
 442   //     for (Int_t ipart=offset; ipart<nMCTracks; ipart++) {
 443
 444   //    AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 445
 446   //    // We don't care about neutrals and non-physical primaries
 447   //    if(mcPart->Charge() == 0) continue;
 448
 449   // PHYSICAL PRIMARY
 450   //    // Get MC vertex
 451     //  TArrayF   vertex;
 452   //    fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
 453   //    Float_t zv = vertex[2];
 454   //    //      Float_t zv = vtxESD->GetZ();
 455   //    // Fill generated histo
 456   //    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
 457
 458   //     }
 459   //   }
 460   // }
 461
 462
 463
 464
 465 }
 466
 467 void   AliAnalysisTaskTriggerStudy::Terminate(Option_t *){
 468   // terminate
 469   // Save output in a more friendly format
 470   fHistoList = dynamic_cast<AliHistoListWrapper*> (GetOutputData(1));
 471   if (!fHistoList){
 472     Printf("ERROR: fHistoList not available");
 473     return;
 474   }
 475   TFile * f = new TFile("trigger_study.root", "recreate");
 476   fHistoList->GetList()->Write();
 477   f->Close();
 478
 479 }
 480
 481 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracklets(const char * name, const char * title){
 482   // Book histo of events vs ntracklets, if needed
 483
 484   TString hname = "hTracklets_";
 485   hname+=name;
 486   hname+=fHistoSuffix;
 487   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 488
 489   if(!h) {
 490     AliInfo(Form("Booking histo %s",hname.Data()));
 491     Bool_t oldStatus = TH1::AddDirectoryStatus();
 492     TH1::AddDirectory(kFALSE);
 493     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 494     h = new TH1F (hname.Data(), title, 4000, -0.5, 4000-0.5);
 495     h->Sumw2();
 496     h->SetXTitle("ntracklets");
 497     fHistoList->GetList()->Add(h);
 498     TH1::AddDirectory(oldStatus);
 499   }
 500   return h;
 501 }
 502
 503
 504 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoSPD1(const char * name, const char * title){
 505   // Book histo of events vs ntracklets, if needed
 506
 507   TString hname = "hSPD1_";
 508   hname+=name;
 509   hname+=fHistoSuffix;
 510   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 511
 512   // const Int_t nbin = 118;
 513   // Float_t bins[nbin] = {0};
 514   // bins[0] = 0;
 515   // for(Int_t ibin = 1; ibin <= nbin; ibin++){
 516   //   if (ibin < 100)
 517   //     bins[ibin] = bins [ibin-1]+1;
 518   //   else if (ibin < 1000)
 519   //     bins[ibin] = bins [ibin-1]+100;
 520   //   else if (ibin < 10000)
 521   //     bins[ibin] = bins [ibin-1]+1000;
 522   // }
 523
 524
 525   if(!h) {
 526     AliInfo(Form("Booking histo %s",hname.Data()));
 527     Bool_t oldStatus = TH1::AddDirectoryStatus();
 528     TH1::AddDirectory(kFALSE);
 529     //    h = new TH1F (hname.Data(), title, nbin, bins);
 530     //h = new TH1F (hname.Data(), title, 100, -0.5, 100-0.5);
 531     h = new TH1F (hname.Data(), title, 10000, -0.5, 10000-0.5);
 532     h->Sumw2();
 533     h->SetXTitle("SPD1");
 534     fHistoList->GetList()->Add(h);
 535     TH1::AddDirectory(oldStatus);
 536   }
 537   return h;
 538 }
 539
 540 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoV0M(const char * name, const char * title){
 541   // Book histo of events vs ntracklets, if needed
 542
 543   TString hname = "hV0M_";
 544   hname+=name;
 545   hname+=fHistoSuffix;
 546   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 547
 548   if(!h) {
 549     AliInfo(Form("Booking histo %s",hname.Data()));
 550     Bool_t oldStatus = TH1::AddDirectoryStatus();
 551     TH1::AddDirectory(kFALSE);
 552     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 553     h = new TH1F (hname.Data(), title, 300, -0.5, 300-0.5);
 554     h->Sumw2();
 555     h->SetXTitle("V0M");
 556     fHistoList->GetList()->Add(h);
 557     TH1::AddDirectory(oldStatus);
 558   }
 559   return h;
 560 }
 561
 562
 563 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracks(const char * name, const char * title){
 564   // Book histo of events vs ntracklets, if needed
 565
 566   TString hname = "hTPCTracks_";
 567   hname+=name;
 568   hname+=fHistoSuffix;
 569   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 570
 571   if(!h) {
 572     AliInfo(Form("Booking histo %s",hname.Data()));
 573     Bool_t oldStatus = TH1::AddDirectoryStatus();
 574     TH1::AddDirectory(kFALSE);
 575     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 576     h = new TH1F (hname.Data(), title, 50, -0.5, 50-0.5);
 577     h->Sumw2();
 578     h->SetXTitle("ntracks");
 579     fHistoList->GetList()->Add(h);
 580     TH1::AddDirectory(oldStatus);
 581   }
 582   return h;
 583 }
 584
 585
 586 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoEta(const char * name, const char * title){
 587   // Book histo of events vs ntracklets, if needed
 588
 589   TString hname = "hEta_";
 590   hname+=name;
 591   hname+=fHistoSuffix;
 592   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 593
 594   if(!h) {
 595     AliInfo(Form("Booking histo %s",hname.Data()));
 596     Bool_t oldStatus = TH1::AddDirectoryStatus();
 597     TH1::AddDirectory(kFALSE);
 598     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 599     h = new TH1F (hname.Data(), title, 20, -2, 2);
 600     h->Sumw2();
 601     h->SetXTitle("#eta");
 602     fHistoList->GetList()->Add(h);
 603     TH1::AddDirectory(oldStatus);
 604   }
 605   return h;
 606 }
 607
 608 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoPt(const char * name, const char * title){
 609   // Book histo of pt distribution, if needed
 610
 611   TString hname = "hPt_";
 612   hname+=name;
 613   hname+=fHistoSuffix;
 614   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 615
 616   if(!h) {
 617     AliInfo(Form("Booking histo %s",hname.Data()));
 618     Bool_t oldStatus = TH1::AddDirectoryStatus();
 619     TH1::AddDirectory(kFALSE);
 620     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 621     h = new TH1F (hname.Data(), title, 100, -0.5, 1-0.5);
 622     h->Sumw2();
 623     h->SetXTitle("p_{T}");
 624     fHistoList->GetList()->Add(h);
 625     TH1::AddDirectory(oldStatus);
 626   }
 627   return h;
 628 }
 629
 630 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationSPDTPCVz(const char * name, const char * title){
 631   // Book histo of events vs ntracklets, if needed
 632
 633   TString hname = "hTPCvsSPD_";
 634   hname+=name;
 635   hname+=fHistoSuffix;
 636   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 637
 638   if(!h) {
 639     AliInfo(Form("Booking histo %s",hname.Data()));
 640     Bool_t oldStatus = TH1::AddDirectoryStatus();
 641     TH1::AddDirectory(kFALSE);
 642     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 643     h = new TH2F (hname.Data(), title, 80, -20, 20,  80, -20, 20);
 644     h->Sumw2();
 645     h->SetXTitle("SPD Vz");
 646     h->SetYTitle("TPC Vz");
 647     fHistoList->GetList()->Add(h);
 648     TH1::AddDirectory(oldStatus);
 649   }
 650   return h;
 651 }
 652
 653 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationContrTPCSPDCls(const char * name, const char * title){
 654   // Book histo of events vs ntracklets, if needed
 655
 656   TString hname = "hContrTPCvsSPDCls_";
 657   hname+=name;
 658   hname+=fHistoSuffix;
 659   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 660
 661   if(!h) {
 662     AliInfo(Form("Booking histo %s",hname.Data()));
 663     Bool_t oldStatus = TH1::AddDirectoryStatus();
 664     TH1::AddDirectory(kFALSE);
 665     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 666     h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
 667     h->Sumw2();
 668     h->SetXTitle("SPD Clusters");
 669     h->SetYTitle("TPC V Contributors");
 670     fHistoList->GetList()->Add(h);
 671     TH1::AddDirectory(oldStatus);
 672   }
 673   return h;
 674 }
 675 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationTrackletsSPDCls(const char * name, const char * title){
 676   // Book histo of events vs ntracklets, if needed
 677
 678   TString hname = "hTrackletsvsSPDCls_";
 679   hname+=name;
 680   hname+=fHistoSuffix;
 681   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 682
 683   if(!h) {
 684     AliInfo(Form("Booking histo %s",hname.Data()));
 685     Bool_t oldStatus = TH1::AddDirectoryStatus();
 686     TH1::AddDirectory(kFALSE);
 687     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
 688     h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
 689     h->Sumw2();
 690     h->SetXTitle("SPD Clusters");
 691     h->SetYTitle("N tracklets");
 692     fHistoList->GetList()->Add(h);
 693     TH1::AddDirectory(oldStatus);
 694   }
 695   return h;
 696 }
 697
 698
 699 void AliAnalysisTaskTriggerStudy::FillTriggerOverlaps (const char * name, const char * title, Bool_t * vdArray){
 700   //Fills a histo with the different trigger statistics in a venn like diagramm. Books it if needed.
 701
 702   // Get or book histo
 703   TString hname = "hTrigStat_";
 704   hname+=name;
 705   hname+=fHistoSuffix;
 706   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
 707
 708   if(!h) {
 709     AliInfo(Form("Booking histo %s",hname.Data()));
 710     Bool_t oldStatus = TH1::AddDirectoryStatus();
 711     TH1::AddDirectory(kFALSE);
 712     Int_t nbins = 0;
 713     for(Int_t ientry = 0; ientry < kNVDEntries; ientry++){
 714       nbins = nbins | (1<<ientry);
 715     }
 716     //    cout << "NBINS " << nbins << endl;
 717     nbins = nbins+1;
 718     h = new TH1I (hname, title, nbins, -0.5, nbins-0.5);
 719     fHistoList->GetList()->Add(h);
 720     TH1::AddDirectory(oldStatus);
 721
 722     // we look at the combinations of n triggers
 723     // We set a bit for each trigger to fill the diagram
 724     // This is much simpler and faster than any recursive function
 725     h->GetXaxis()->SetBinLabel(1,"NONE");
 726     for(Int_t ibin = 1; ibin < nbins; ibin++){
 727       TString binname = "";
 728       Bool_t first = kTRUE;
 729       for(Int_t ivdentry = 0; ivdentry < kNVDEntries; ivdentry++){
 730         if (ibin & (1<<ivdentry)) {
 731           if(!first) binname += " & ";
 732           binname += kVDNames[ivdentry];
 733           first=kFALSE;
 734         }
 735       }
 736       h->GetXaxis()->SetBinLabel(ibin+1,binname.Data());
 737     }
 738
 739   }
 740
 741   UInt_t mask = 0;
 742   for(Int_t ivdentry = 0; ivdentry < kNVDEntries; ivdentry++){
 743     if(vdArray[ivdentry]) {
 744       mask  = mask | (1<<ivdentry);
 745       //      cout << " 1 "   ;
 746     } //else cout << " 0 ";
 747   }
 748   //  cout << hex << " = " << mask << endl;
 749
 750   h->Fill(mask);
 751
 752 }