]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/PHOSTasks/PHOS_PbPb/AliAnalysisTaskPi0Flow.cxx
added FillHistogram for 3d Hists with weight
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / AliAnalysisTaskPi0Flow.cxx
index 9b624d2813c23bc15bd9a699321ddf2ca07f0f8f..05e94d100cd8026da3ab8bdaafc00ff82fd99de5 100644 (file)
@@ -88,7 +88,6 @@ AliAnalysisTaskPi0Flow::AliAnalysisTaskPi0Flow(const char *name, Period period)
   fEvent(0x0),
   fEventESD(0x0),
   fEventAOD(0x0),
-  fMCStack(0x0),
   fRunNumber(-999),
   fInternalRunNumber(0),
   fPHOSGeo(0),
@@ -119,6 +118,9 @@ AliAnalysisTaskPi0Flow::AliAnalysisTaskPi0Flow(const char *name, Period period)
   TArrayI centNMixed(nbins, nMixed);
   SetCentralityBinning(centEdges, centNMixed);
   
+  for(int mod=1; mod <= kNMod; ++mod)
+    fModuleEnabled[mod-1] = kTRUE;
+
   for(Int_t i=0;i<kNCenBins;i++){
     for(Int_t j=0;j<2; j++)
       for(Int_t k=0; k<2; k++) {
@@ -339,6 +341,8 @@ void AliAnalysisTaskPi0Flow::UserCreateOutputObjects()
     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
   }
   
+  this->MakeMCHistograms();
+  
   // Setup photon lists
   Int_t kapacity = kNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
   fCaloPhotonsPHOSLists = new TObjArray(kapacity);
@@ -347,6 +351,17 @@ void AliAnalysisTaskPi0Flow::UserCreateOutputObjects()
   PostData(1, fOutputContainer);
 }
 
+void AliAnalysisTaskPi0Flow::MakeMCHistograms()
+{
+  //empty, MC extensions should override
+}
+
+void AliAnalysisTaskPi0Flow::DoMC()
+{
+  //empty, MC extensions should override
+}
+
+
 //________________________________________________________________________
 void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
 {
@@ -358,7 +373,6 @@ void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
   fEvent = GetEvent();
   fEventESD = dynamic_cast<AliESDEvent*> (fEvent);
   fEventAOD = dynamic_cast<AliAODEvent*> (fEvent);
-  fMCStack = GetMCStack();
   LogProgress(0);
 
 
@@ -447,7 +461,10 @@ void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
   ConsiderPi0sMix();
   LogProgress(10);
   
-  // Step 11: Update lists
+  // Step 11: MC
+  this->DoMC();
+
+  // Step 12: Update lists
   UpdateLists();
   LogProgress(11);
 
@@ -479,6 +496,14 @@ void AliAnalysisTaskPi0Flow::SetCentralityBinning(const TArrayD& edges, const TA
   fCentNMixed = nMixed;
 }
 
+//_____________________________________________________________________________
+void AliAnalysisTaskPi0Flow::SetEnablePHOSModule(int module, Bool_t enable)
+{
+  if( module < 1 || kNMod < module )
+    AliFatal(Form("PHOS Module must be between 1 and %i", kNMod));
+  else
+    fModuleEnabled[module-1] = enable;
+}
 
 
 //________________________________________________________________________
@@ -599,6 +624,8 @@ void AliAnalysisTaskPi0Flow::SelectPhotonClusters()
     Int_t mod  = relId[0] ;
     Int_t cellX = relId[2];
     Int_t cellZ = relId[3] ;
+    if ( ! fModuleEnabled[mod-1] )
+      continue;
     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
       continue ; // reject if not.
 
@@ -663,6 +690,7 @@ void AliAnalysisTaskPi0Flow::SelectPhotonClusters()
 
     fCaloPhotonsPHOS->Add(new  AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E()) );
     AliCaloPhoton * ph = (AliCaloPhoton*) fCaloPhotonsPHOS->At( fCaloPhotonsPHOS->GetLast() );
+    ph->SetCluster(clu);
 
     ph->SetModule(mod) ;
     lorentzMomentum*= ecore/lorentzMomentum.E() ;
@@ -672,7 +700,8 @@ void AliAnalysisTaskPi0Flow::SelectPhotonClusters()
     //Evaluate CoreDispersion
     Double_t m02=0.,m20=0. ;
     EvalCoreLambdas(clu, cells, m02, m20) ;
-    ph->SetDisp2Bit(TestCoreLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
+    ph->SetDisp2Bit(TestCoreLambda(clu->E(),m20,m02)) ; //Correct order m20,m02
+//    ph->SetDisp2Bit(TestCoreLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
     if(ph->IsDispOK()){
       FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.);
     }
@@ -1417,6 +1446,20 @@ void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t
   AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
 }
 
+void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z, Double_t w) const{
+  //Fills 1D histograms with key
+  TObject * obj = fOutputContainer->FindObject(key);
+  
+  TH3 * th3 = dynamic_cast<TH3*> (obj);
+  if(th3) {
+    th3->Fill(x, y, z, w) ;
+    return;
+  }
+  
+  AliError(Form("can not find histogram (of instance TH3) <%s> ",key)) ;
+}
+
+
 //_____________________________________________________________________________
 AliVEvent* AliAnalysisTaskPi0Flow::GetEvent()
 {
@@ -1430,17 +1473,17 @@ AliVEvent* AliAnalysisTaskPi0Flow::GetEvent()
 
 
 //___________________________________________________________________________
-AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
-{
-  fMCStack = 0;
-  AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
-  if(eventHandler){
-    AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
-    if( mcEventHandler)
-      fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
-  }
-  return fMCStack;
-}
+// AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
+// {
+//   fMCStack = 0;
+//   AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
+//   if(eventHandler){
+//     AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
+//     if( mcEventHandler)
+//       fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
+//   }
+//   return fMCStack;
+// }
 
 //___________________________________________________________________________
 Int_t AliAnalysisTaskPi0Flow::GetCentralityBin(Float_t centralityV0M)
@@ -1472,7 +1515,8 @@ Int_t AliAnalysisTaskPi0Flow::GetRPBin()
 {
   Double_t averageRP;
   if(fHaveTPCRP)
-    averageRP = (fRPV0A+fRPV0C+fRP) /3.;
+    averageRP = fRP ;// If possible, it is better to have EP bin from TPC
+                     // to have similar events for miximng (including jets etc)   (fRPV0A+fRPV0C+fRP) /3.;
   else
     averageRP = (fRPV0A+fRPV0C) /2.;
 
@@ -1841,6 +1885,181 @@ Int_t AliAnalysisTaskPi0Flow::ConvertToInternalRunNumber(Int_t run){
     default : return 199;
     }
   }
+  if( kLHC13 == fPeriod ) {
+    switch(run){
+    case  195344 : return 1;
+    case  195346 : return 2;
+    case  195351 : return 3;
+    case  195389 : return 4;
+    case  195390 : return 5;
+    case  195391 : return 6;
+    case  195478 : return 7;
+    case  195479 : return 8;
+    case  195480 : return 9;
+    case  195481 : return 10;
+    case  195482 : return 11;
+    case  195483 : return 12;
+    case  195529 : return 13;
+    case  195531 : return 14;
+    case  195532 : return 15;
+    case  195566 : return 16;
+    case  195567 : return 17;
+    case  195568 : return 18;
+    case  195592 : return 19;
+    case  195593 : return 20;
+    case  195596 : return 21;
+    case  195633 : return 22;
+    case  195635 : return 23;
+    case  195644 : return 24;
+    case  195673 : return 25;
+    case  195675 : return 26;
+    case  195676 : return 27;
+    case  195677 : return 28;
+    case  195681 : return 29;
+    case  195682 : return 30;
+    case  195720 : return 31;
+    case  195721 : return 32;
+    case  195722 : return 33;
+    case  195724 : return 34;
+    case  195725 : return 34;
+    case  195726 : return 35;
+    case  195727 : return 36;
+    case  195760 : return 37;
+    case  195761 : return 38;
+    case  195765 : return 39;
+    case  195767 : return 40;
+    case  195783 : return 41;
+    case  195787 : return 42;
+    case  195826 : return 43;
+    case  195827 : return 44;
+    case  195829 : return 45;
+    case  195830 : return 46;
+    case  195831 : return 47;
+    case  195867 : return 48;
+    case  195869 : return 49;
+    case  195871 : return 50;
+    case  195872 : return 51;
+    case  195873 : return 52;
+    case  195935 : return 53;
+    case  195949 : return 54;
+    case  195950 : return 55;
+    case  195954 : return 56;
+    case  195955 : return 57;
+    case  195958 : return 58;
+    case  195989 : return 59;
+    case  195994 : return 60;
+    case  195998 : return 61;
+    case  196000 : return 62;
+    case  196006 : return 63;
+    case  196085 : return 64;
+    case  196089 : return 65;
+    case  196090 : return 66;
+    case  196091 : return 67;
+    case  196099 : return 68;
+    case  196105 : return 69;
+    case  196107 : return 70;
+    case  196185 : return 71;
+    case  196187 : return 72;
+    case  196194 : return 73;
+    case  196197 : return 74;
+    case  196199 : return 75;
+    case  196200 : return 76;
+    case  196201 : return 77;
+    case  196203 : return 78;
+    case  196208 : return 79;
+    case  196214 : return 80;
+    case  196308 : return 81;
+    case  196309 : return 82;
+    case  196310 : return 83;
+    case  196311 : return 84;
+    case  196433 : return 85;
+    case  196474 : return 86;
+    case  196475 : return 87;
+    case  196477 : return 88;
+    case  196528 : return 89;
+    case  196533 : return 90;
+    case  196535 : return 91;
+    case  196563 : return 92;
+    case  196564 : return 93;
+    case  196566 : return 94;
+    case  196568 : return 95;
+    case  196601 : return 96;
+    case  196605 : return 97;
+    case  196608 : return 98;
+    case  196646 : return 99;
+    case  196648 : return 100;
+    case  196701 : return 101;
+    case  196702 : return 102;
+    case  196703 : return 103;
+    case  196706 : return 104;
+    case  196714 : return 105;
+    case  196720 : return 106;
+    case  196721 : return 107;
+    case  196722 : return 108;
+    case  196772 : return 109;
+    case  196773 : return 110;
+    case  196774 : return 111;
+    case  196869 : return 112;
+    case  196870 : return 113;
+    case  196874 : return 114;
+    case  196876 : return 115;
+    case  196965 : return 116;
+    case  196967 : return 117;
+    case  196972 : return 118;
+    case  196973 : return 119;
+    case  196974 : return 120;
+    case  197003 : return 121;
+    case  197011 : return 122;
+    case  197012 : return 123;
+    case  197015 : return 124;
+    case  197027 : return 125;
+    case  197031 : return 126;
+    case  197089 : return 127;
+    case  197090 : return 128;
+    case  197091 : return 129;
+    case  197092 : return 130;
+    case  197094 : return 131;
+    case  197098 : return 132;
+    case  197099 : return 133;
+    case  197138 : return 134;
+    case  197139 : return 135;
+    case  197142 : return 136;
+    case  197143 : return 137;
+    case  197144 : return 138;
+    case  197145 : return 139;
+    case  197146 : return 140;
+    case  197147 : return 140;
+    case  197148 : return 141;
+    case  197149 : return 142;
+    case  197150 : return 143;
+    case  197152 : return 144;
+    case  197153 : return 145;
+    case  197184 : return 146;
+    case  197189 : return 147;
+    case  197247 : return 148;
+    case  197248 : return 149;
+    case  197254 : return 150;
+    case  197255 : return 151;
+    case  197256 : return 152;
+    case  197258 : return 153;
+    case  197260 : return 154;
+    case  197296 : return 155;
+    case  197297 : return 156;
+    case  197298 : return 157;
+    case  197299 : return 158;
+    case  197300 : return 159;
+    case  197302 : return 160;
+    case  197341 : return 161;
+    case  197342 : return 162;
+    case  197348 : return 163;
+    case  197349 : return 164;
+    case  197351 : return 165;
+    case  197386 : return 166;
+    case  197387 : return 167;
+    case  197388 : return 168;
+    default : return 199;
+    }
+  }
   if((fPeriod == kUndefinedPeriod) && (fDebug >= 1) ) {
     AliWarning("Period not defined");
   }
@@ -2548,7 +2767,7 @@ void  AliAnalysisTaskPi0Flow::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells *
   const Int_t mulDigit=clu->GetNCells() ;
   Double_t xc[mulDigit] ;
   Double_t zc[mulDigit] ;
-  Double_t ei[mulDigit] ;
+  Double_t wi[mulDigit] ;
   Double_t x = 0 ;
   Double_t z = 0 ;
   const Double_t logWeight=4.5 ;
@@ -2561,9 +2780,11 @@ void  AliAnalysisTaskPi0Flow::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells *
     fPHOSGeo->RelPosInModule(relid, xi, zi);
     xc[iDigit]=xi ;
     zc[iDigit]=zi ;
-    ei[iDigit] = elist[iDigit]*cells->GetCellAmplitude(absId) ;
-    if (clu->E()>0 && ei[iDigit]>0) {
-      Float_t w = TMath::Max( 0., logWeight + TMath::Log( ei[iDigit] / clu->E() ) ) ;
+    Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
+    wi[iDigit]=0. ;
+    if (clu->E()>0 && ei>0) {
+      wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
+      Double_t w=wi[iDigit];
       x    += xc[iDigit] * w ;
       z    += zc[iDigit] * w ;
       wtot += w ;
@@ -2581,8 +2802,8 @@ void  AliAnalysisTaskPi0Flow::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells *
   Double_t xCut = 0. ;
   Double_t zCut = 0. ;
   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
-    if (clu->E()>0 && ei[iDigit]>0.) {
-        Double_t w = TMath::Max( 0., logWeight + TMath::Log( ei[iDigit] / clu->E() ) ) ;
+    Double_t w=wi[iDigit];
+    if (w>0.) {
         Double_t xi= xc[iDigit] ;
         Double_t zi= zc[iDigit] ;
        if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){