PWGGA/PHOSTasks/PHOSCalib/AliAnalysisTaskPi0Calib.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TObjArray.h"
   4 #include "TF1.h"
   5 #include "TFile.h"
   6 #include "TH1F.h"
   7 #include "TH2F.h"
   8 #include "TH2I.h"
   9 #include "TH3F.h"
  10 #include "TParticle.h"
  11 #include "TCanvas.h"
  12 #include "TStyle.h"
  13 #include "TRandom.h"
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliAnalysisTaskSE.h"
  17 #include "AliAnalysisTaskPi0Calib.h"
  18 #include "AliCaloPhotonC.h"
  19 #include "AliPHOSGeometry.h"
  20 #include "AliPHOSEsdCluster.h"
  21 #include "AliPHOSCalibData.h"
  22 #include "AliESDEvent.h"
  23 #include "AliESDCaloCells.h"
  24 #include "AliESDVertex.h"
  25 #include "AliLog.h"
  26 #include "AliPID.h"
  27 #include "AliCDBManager.h"
  28
  29 // Analysis task to fill histograms with PHOS ESD clusters and cells
  30 // Authors: Dmitri Peressounko
  31 // Date   : 28.05.2011
  32
  33 ClassImp(AliAnalysisTaskPi0Calib)
  34
  35 //________________________________________________________________________
  36 Double_t rnlin(Double_t *x, Double_t *par)
  37 {
  38   //a = par[0], b = par[1].
  39   //1+a*exp(-e/b)
  40
  41  return 0.0241+1.0504*x[0]+0.000249*x[0]*x[0] ;
  42
  43   Double_t gcbCorr = 0.0241+1.0504*x[0]+0.000249*x[0]*x[0] ;
  44 //  Double_t bvpCorr = 1+par[0]*TMath::Exp(-x[0]/par[1]);
  45 //  Double_t bvpCorr = 1.005*(1+par[0]*par[1]*par[1]/(x[0]*x[0]+par[1]*par[1])) ;
  46   Double_t bvpCorr = 1.007*(1+par[0]*par[1]*par[1]/(x[0]*x[0]+par[1]*par[1])) ;
  47   return gcbCorr*bvpCorr;
  48
  49 }
  50 //________________________________________________________________________
  51 Double_t nonlin(Double_t e){
  52  return 0.0241+1.0504*e+0.000249*e*e ;
  53 }
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskPi0Calib::AliAnalysisTaskPi0Calib(const char *name)
  57   : AliAnalysisTaskSE(name),
  58     fOutputContainer(0),
  59     fPHOSEvents(0),
  60     fPHOSEvent(0),
  61     fPHOSCalibData(0),
  62     fPHOSGeo(0),
  63     fEventCounter(0)
  64 {
  65   // Constructor
  66
  67   // Output slots #0 write into a TH1 container
  68   DefineOutput(1,TList::Class());
  69
  70   // Set bad channel map
  71   char key[55] ;
  72   for(Int_t i=0; i<6; i++){
  73     sprintf(key,"PHOS_BadMap_mod%d",i) ;
  74     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
  75     fPHOSAmp[i]=0x0 ;
  76   }
  77   // Initialize the PHOS geometry
  78   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
  79
  80   //We have to apply re-calibration for pass1 LCH10h
  81   // Initialize decalibration factors in the form of the OCDB object
  82
  83
  84   AliCDBManager * man = AliCDBManager::Instance();
  85   man->SetRun(140000) ;
  86   man->SetDefaultStorage("local://OCDB");
  87   fPHOSCalibData = new AliPHOSCalibData();
  88
  89
  90   // Initialize non-linrarity correction
  91 //  fNonLinCorr = new TF1("nonlib",rnlin,0.,40.,0);
  92
  93   fOutputContainer = new TList();
  94   fOutputContainer->SetOwner(kTRUE);
  95
  96 }
  97 //________________________________________________________________________
  98 AliAnalysisTaskPi0Calib::~AliAnalysisTaskPi0Calib(){
  99
 100   for(Int_t i=0; i<6; i++){
 101     if(fPHOSBadMap[i]){
 102       delete fPHOSBadMap[i] ;
 103       fPHOSBadMap[i]=0x0 ;
 104     }
 105     if(fPHOSAmp[i]){
 106       delete fPHOSAmp[i] ;
 107       fPHOSAmp[i]=0x0 ;
 108     }
 109   }
 110 }
 111
 112 //________________________________________________________________________
 113 void AliAnalysisTaskPi0Calib::UserCreateOutputObjects()
 114 {
 115   // Create histograms
 116
 117 //  AliCDBManager * man = AliCDBManager::Instance();
 118 //  man->SetRun(140000) ;
 119 //  man->SetDefaultStorage("local://OCDB");
 120 //  fPHOSCalibData = new AliPHOSCalibData();
 121
 122 //  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
 123
 124
 125
 126   // Called once
 127   const Int_t nRuns=200 ;
 128
 129   // ESD histograms
 130 //  if(fOutputContainer != NULL){
 131 //    delete fOutputContainer;
 132 //  }
 133 //  fOutputContainer = new TList();
 134 //  fOutputContainer->SetOwner(kTRUE);
 135
 136   //========QA histograms=======
 137
 138   //Event selection
 139   fOutputContainer->Add(new TH1F("hSelEvents","Event selection", 10,0.,10.)) ;
 140
 141   //vertex distribution
 142   fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
 143
 144   fOutputContainer->Add(new TH3F("hMggAllM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 145   fOutputContainer->Add(new TH3F("hMggAllM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 146   fOutputContainer->Add(new TH3F("hMggAllM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 147
 148   fOutputContainer->Add(new TH3F("hMggDispM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 149   fOutputContainer->Add(new TH3F("hMggDispM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 150   fOutputContainer->Add(new TH3F("hMggDispM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 151
 152   fOutputContainer->Add(new TH3F("hMggBothM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 153   fOutputContainer->Add(new TH3F("hMggBothM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 154   fOutputContainer->Add(new TH3F("hMggBothM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 155
 156   fOutputContainer->Add(new TH3F("hMiMggAllM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 157   fOutputContainer->Add(new TH3F("hMiMggAllM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 158   fOutputContainer->Add(new TH3F("hMiMggAllM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 159
 160   fOutputContainer->Add(new TH3F("hMiMggDispM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 161   fOutputContainer->Add(new TH3F("hMiMggDispM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 162   fOutputContainer->Add(new TH3F("hMiMggDispM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 163
 164   fOutputContainer->Add(new TH3F("hMggOldM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 165   fOutputContainer->Add(new TH3F("hMggOldM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 166   fOutputContainer->Add(new TH3F("hMggOldM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 167
 168   fOutputContainer->Add(new TH3F("hMggTOFM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 169   fOutputContainer->Add(new TH3F("hMggTOFM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 170   fOutputContainer->Add(new TH3F("hMggTOFM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,0.,0.5));
 171
 172
 173   Int_t nPtPhot = 300 ;
 174   Double_t ptPhotMax = 30 ;
 175   Int_t nM       = 500;
 176   Double_t mMin  = 0.0;
 177   Double_t mMax  = 1.0;
 178   fOutputContainer->Add(new TH2F("hPi0M11","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 179   fOutputContainer->Add(new TH2F("hPi0M22","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 180   fOutputContainer->Add(new TH2F("hPi0M33","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 181   fOutputContainer->Add(new TH2F("hPi0M12","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 182   fOutputContainer->Add(new TH2F("hPi0M23","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 183   fOutputContainer->Add(new TH2F("hPi0M13","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 184
 185   fOutputContainer->Add(new TH2F("hMiPi0M11","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 186   fOutputContainer->Add(new TH2F("hMiPi0M22","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 187   fOutputContainer->Add(new TH2F("hMiPi0M33","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 188   fOutputContainer->Add(new TH2F("hMiPi0M12","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 189   fOutputContainer->Add(new TH2F("hMiPi0M23","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 190   fOutputContainer->Add(new TH2F("hMiPi0M13","All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 191
 192   fOutputContainer->Add(new TH3F("hTOFM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,-10.e-6,10.e-6));
 193   fOutputContainer->Add(new TH3F("hTOFM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,-10.e-6,10.e-6));
 194   fOutputContainer->Add(new TH3F("hTOFM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,-10.e-6,10.e-6));
 195   fOutputContainer->Add(new TH3F("hTOFFineM1","M1" ,64,0.5,64.5, 56,0.5,56.5,200,-100.e-9,100.e-9));
 196   fOutputContainer->Add(new TH3F("hTOFFineM2","M2" ,64,0.5,64.5, 56,0.5,56.5,200,-100.e-9,100.e-9));
 197   fOutputContainer->Add(new TH3F("hTOFFineM3","M3" ,64,0.5,64.5, 56,0.5,56.5,200,-100.e-9,100.e-9));
 198
 199   fOutputContainer->Add(new TH2F("hTOFvsEM1","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 200   fOutputContainer->Add(new TH2F("hTOFvsEM2","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 201   fOutputContainer->Add(new TH2F("hTOFvsEM3","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 202
 203   fOutputContainer->Add(new TH2F("hTOFAvvsEM1","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 204   fOutputContainer->Add(new TH2F("hTOFAvvsEM2","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 205   fOutputContainer->Add(new TH2F("hTOFAvvsEM3","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 206
 207   fOutputContainer->Add(new TH2F("hTOFvsECluM1","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 208   fOutputContainer->Add(new TH2F("hTOFvsECluM2","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 209   fOutputContainer->Add(new TH2F("hTOFvsECluM3","All clusters",200,-100.e-9,100.e-9,100,0.,20.));
 210
 211
 212   PostData(1, fOutputContainer);
 213
 214 }
 215
 216 //________________________________________________________________________
 217 void AliAnalysisTaskPi0Calib::UserExec(Option_t *)
 218 {
 219   // Main loop, called for each event
 220   // Analyze ESD/AOD
 221
 222   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 223   if (!event) {
 224     Printf("ERROR: Could not retrieve event");
 225     PostData(1, fOutputContainer);
 226     return;
 227   }
 228   FillHistogram("hSelEvents",0.5) ;
 229
 230   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 231 //  char skey[55] ;
 232   if(fEventCounter == 0) {
 233     for(Int_t mod=0; mod<5; mod++) {
 234       if(!event->GetPHOSMatrix(mod)) continue;
 235       fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
 236       Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
 237     }
 238     fEventCounter++ ;
 239   }
 240
 241   if(fPHOSCalibData->GetADCchannelEmc(3,15,26)==1.){
 242     for(Int_t module=1; module<4; module++) {
 243       for(Int_t ix=1; ix<=64; ix++) {
 244         for(Int_t iz=1; iz<=56; iz++) {
 245
 246           Float_t corr = fPHOSAmp[module]->GetBinContent(ix,iz) ;
 247             fPHOSCalibData->SetADCchannelEmc(module,iz,ix,corr);
 248         }
 249       }
 250     }
 251   }
 252
 253
 254   // Checks if we have a primary vertex
 255   // Get primary vertices form ESD
 256   const AliESDVertex *esdVertex5 = event->GetPrimaryVertex();
 257
 258   Double_t vtx0[3] = {0,0,0}; // don't rely on ESD vertex, assume (0,0,0)
 259   Double_t vtx5[3] ={0.,0.,0.};
 260
 261   vtx5[0] = esdVertex5->GetX();
 262   vtx5[1] = esdVertex5->GetY();
 263   vtx5[2] = esdVertex5->GetZ();
 264
 265
 266   FillHistogram("hZvertex",esdVertex5->GetZ());
 267 /*
 268   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 269     PostData(1, fOutputContainer);
 270     return;
 271   }
 272 */
 273   FillHistogram("hSelEvents",1.5) ;
 274
 275   if(!fPHOSEvents)
 276     fPHOSEvents=new TList() ;
 277   TList * prevPHOS = fPHOSEvents ;
 278
 279   if(fPHOSEvent)
 280     fPHOSEvent->Clear() ;
 281   else
 282     fPHOSEvent = new TClonesArray("AliCaloPhotonC",200) ;
 283
 284   //For re-calibration
 285   const Double_t logWeight=4.5 ;
 286
 287   TVector3 vertex(vtx5);
 288
 289   Int_t multClust = event->GetNumberOfCaloClusters();
 290   Int_t inPHOS=0; //,inMod1=0, inMod2=0, inMod3=0 ;
 291 //   Double_t avrgEm1=0,avrgEm2=0,avrgEm3=0; //average cluster energy
 292
 293   AliESDCaloCells * cells = event->GetPHOSCells() ;
 294 //   char key[55];
 295
 296   TVector3 globalPos;
 297   TVector3 localPos ;
 298   TLorentzVector p12 ;
 299   for (Int_t i=0; i<multClust; i++) {
 300     AliESDCaloCluster *clu = event->GetCaloCluster(i);
 301     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 302
 303     Float_t  position[3];
 304     clu->GetPosition(position);
 305     TVector3 global(position) ;
 306     Int_t relId[4] ;
 307     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 308     Int_t mod  = relId[0] ;
 309     Int_t cellX = relId[2];
 310     Int_t cellZ = relId[3] ;
 311 //    if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 312 //      continue ;
 313
 314
 315     //Apply re-Calibreation
 316     AliPHOSEsdCluster cluPHOS1(*clu);
 317     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
 318     cluPHOS1.EvalAll(logWeight,vertex);         // recalculate the cluster parameters
 319     cluPHOS1.SetE(nonlin(cluPHOS1.E()));// Users's nonlinearity
 320     if(cluPHOS1.E()<0.3) continue;
 321     if(clu->GetNCells()<3) continue ;
 322     if(clu->GetM02()<0.2)   continue ;
 323
 324     //TOF recalibration
 325     Double_t dt=0.;
 326     Double_t tof = clu->GetTOF()-dt ;
 327
 328     FillHistogram(Form("hTOFM%d",mod),float(cellX),float(cellZ),tof);
 329     if(clu->E()>1.){
 330       FillHistogram(Form("hTOFFineM%d",mod),float(cellX),float(cellZ),tof);
 331     }
 332     FillHistogram(Form("hTOFvsEM%d",mod),tof,clu->E());
 333
 334 /*
 335     //Calculate TOF within cluster
 336     Double_t tofAv = 0. ;
 337     Double_t wtot=0. ;
 338     for(Int_t iDigit=0; iDigit<clu->GetNCells(); iDigit++) {
 339       Int_t absId=clu->GetCellAbsId(iDigit);
 340       Double_t edig = cells->GetCellAmplitude(absId) ;
 341       Double_t tau = cells->GetCellTime(absId) ;
 342       Int_t relid[4] ;
 343       fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
 344       if(fPHOSTOF[relid[0]]){
 345         tau -= fPHOSTOF[relid[0]]->GetBinContent(relid[2],relid[3]) ;
 346       }
 347       sprintf(key,"hTOFvsECluM%d",mod) ;
 348       FillHistogram(key,tau-tof,edig);
 349
 350       if(edig>0.4 && TMath::Abs(tau)<10.e-9){
 351         Double_t wi=2.73+11.8*TMath::Exp(-edig/0.45) ;
 352         tofAv+=tau/wi/wi ;
 353         wtot+=1./wi/wi ;
 354       }
 355     }
 356     if(wtot>0.)
 357       tofAv/=wtot ;
 358     else
 359       tofAv=tof ;
 360     FillHistogram(Form("hTOFAvvsEM%d",mod),tofAv,clu->E());
 361 */
 362
 363
 364     TLorentzVector pv1 ;
 365     cluPHOS1.GetMomentum(pv1,vtx0);
 366     if(inPHOS>=fPHOSEvent->GetSize()){
 367       fPHOSEvent->Expand(inPHOS+50) ;
 368     }
 369     new((*fPHOSEvent)[inPHOS]) AliCaloPhotonC(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 370     AliCaloPhotonC * ph = (AliCaloPhotonC*)fPHOSEvent->At(inPHOS) ;
 371     ph->SetModule(mod) ;
 372
 373     clu->GetMomentum(pv1 ,vtx0);
 374     ph->SetMomV2(&pv1) ;
 375     ph->SetNCells(clu->GetNCells());
 376     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
 377     ph->SetTOFBit(TestTOF(clu->E(),tof)) ;
 378     //Track matching
 379     Double_t dx=clu->GetTrackDx() ;
 380     Double_t dz=clu->GetTrackDz() ;
 381     Bool_t cpvBit=kTRUE ; //No track matched by default
 382     Bool_t cpvBit2=kTRUE ; //More Strict criterion
 383     TArrayI * itracks = clu->GetTracksMatched() ;
 384     if(itracks->GetSize()>0){
 385       Int_t iTr = itracks->At(0);
 386       if(iTr>=0 && iTr<event->GetNumberOfTracks()){
 387         AliESDtrack *track = event->GetTrack(iTr) ;
 388         Double_t pt = track->Pt() ;
 389         Short_t charge = track->Charge() ;
 390         Double_t r=TestCPV(dx, dz, pt,charge) ;
 391         cpvBit=(r>2.) ;
 392         cpvBit2=(r>4.) ;
 393       }
 394     }
 395     ph->SetCPVBit(cpvBit) ;
 396     ph->SetCPVBit2(cpvBit2) ;
 397     ph->SetEMCx(float(cellX)) ;
 398     ph->SetEMCz(float(cellZ)) ;
 399     ph->SetLambdas(clu->GetM20(),clu->GetM02()) ;
 400     inPHOS++ ;
 401   }
 402
 403
 404   TLorentzVector p12old ;
 405   char key[255] ;
 406   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 407     AliCaloPhotonC * ph1=(AliCaloPhotonC*)fPHOSEvent->At(i1) ;
 408     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 409       AliCaloPhotonC * ph2=(AliCaloPhotonC*)fPHOSEvent->At(i2) ;
 410       p12  = *ph1  + *ph2;
 411       if(p12.Pt()<2.)
 412         continue ;
 413       p12old = *(ph1->GetMomV2()) + *(ph2->GetMomV2()) ;
 414
 415           if(ph1->Module()==1 && ph2->Module()==1)
 416             FillHistogram("hPi0M11",p12.M(),p12.Pt() );
 417           else if(ph1->Module()==2 && ph2->Module()==2)
 418             FillHistogram("hPi0M22",p12.M(),p12.Pt() );
 419           else if(ph1->Module()==3 && ph2->Module()==3)
 420             FillHistogram("hPi0M33",p12.M(),p12.Pt() );
 421           else if(ph1->Module()==1 && ph2->Module()==2)
 422             FillHistogram("hPi0M12",p12.M(),p12.Pt() );
 423           else if(ph1->Module()==1 && ph2->Module()==3)
 424             FillHistogram("hPi0M13",p12.M(),p12.Pt() );
 425           else if(ph1->Module()==2 && ph2->Module()==3)
 426             FillHistogram("hPi0M23",p12.M(),p12.Pt() );
 427
 428       if(ph1->Module()!=ph2->Module())
 429         continue ;
 430
 431       FillHistogram(Form("hMggAllM%d",ph1->Module()),ph1->EMCx(),ph1->EMCz(),p12.M() );
 432       FillHistogram(Form("hMggAllM%d",ph2->Module()),ph2->EMCx(),ph2->EMCz(),p12.M() );
 433
 434       if(ph1->IsTOFOK() && ph2->IsTOFOK()){
 435         sprintf(key,"hMggTOFM%d",ph1->Module()) ;
 436         FillHistogram(Form("hMggTOFM%d",ph1->Module()),ph1->EMCx(),ph1->EMCz(),p12.M() );
 437         FillHistogram(Form("hMggTOFM%d",ph1->Module()),ph2->EMCx(),ph2->EMCz(),p12.M() );
 438       }
 439
 440       if(ph1->IsDispOK() && ph2->IsDispOK()){
 441         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 442           sprintf(key,"hMggDispM%d",ph1->Module()) ;
 443           FillHistogram(key,ph1->EMCx(),ph1->EMCz(),p12.M() );
 444           sprintf(key,"hMggDispM%d",ph2->Module()) ;
 445           FillHistogram(key,ph2->EMCx(),ph2->EMCz(),p12.M() );
 446           if(ph1->IsTOFOK() && ph2->IsTOFOK()){
 447             sprintf(key,"hMggBothM%d",ph1->Module()) ;
 448             FillHistogram(key,ph1->EMCx(),ph1->EMCz(),p12.M() );
 449             sprintf(key,"hMggBothM%d",ph2->Module()) ;
 450             FillHistogram(key,ph2->EMCx(),ph2->EMCz(),p12.M() );
 451           }
 452           sprintf(key,"hMggOldM%d",ph1->Module()) ;
 453           FillHistogram(key,ph1->EMCx(),ph1->EMCz(),p12old.M() );
 454           sprintf(key,"hMggOldM%d",ph2->Module()) ;
 455           FillHistogram(key,ph2->EMCx(),ph2->EMCz(),p12old.M() );
 456         }
 457       }
 458
 459     } // end of loop i2
 460   } // end of loop i1
 461
 462
 463   //now mixed
 464   for (Int_t i1=0; i1<inPHOS; i1++) {
 465     AliCaloPhotonC * ph1=(AliCaloPhotonC*)fPHOSEvent->At(i1) ;
 466     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 467       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 468       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 469         AliCaloPhotonC * ph2=(AliCaloPhotonC*)mixPHOS->At(i2) ;
 470         p12  = *ph1  + *ph2;
 471         if(p12.Pt()<2.)
 472           continue ;
 473
 474       if(ph1->Module()!=ph2->Module())
 475         continue ;
 476
 477             if(ph1->Module()==1 && ph2->Module()==1)
 478               FillHistogram("hMiPi0M11",p12.M(),p12.Pt() );
 479             else if(ph1->Module()==2 && ph2->Module()==2)
 480               FillHistogram("hMiPi0M22",p12.M(),p12.Pt() );
 481             else if(ph1->Module()==3 && ph2->Module()==3)
 482               FillHistogram("hMiPi0M33",p12.M(),p12.Pt() );
 483             else if(ph1->Module()==1 && ph2->Module()==2)
 484               FillHistogram("hMiPi0M12",p12.M(),p12.Pt() );
 485             else if(ph1->Module()==1 && ph2->Module()==3)
 486               FillHistogram("hMiPi0M13",p12.M(),p12.Pt() );
 487             else if(ph1->Module()==2 && ph2->Module()==3)
 488               FillHistogram("hMiPi0M23",p12.M(),p12.Pt() );
 489
 490
 491         sprintf(key,"hMiMggAllM%d",ph1->Module()) ;
 492         FillHistogram(key,ph1->EMCx(),ph1->EMCz(),p12.M() );
 493         FillHistogram(key,ph2->EMCx(),ph2->EMCz(),p12.M() );
 494
 495         if(ph1->IsDispOK() && ph2->IsDispOK()){
 496           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 497             sprintf(key,"hMiMggDispM%d",ph1->Module()) ;
 498             FillHistogram(key,ph1->EMCx(),ph1->EMCz(),p12.M() );
 499             FillHistogram(key,ph2->EMCx(),ph2->EMCz(),p12.M() );
 500           }
 501         }
 502
 503       } // end of loop i2
 504     }
 505   } // end of loop i1
 506
 507
 508   //Now we either add current events to stack or remove
 509   //If no photons in current event - no need to add it to mixed
 510   if(fPHOSEvent->GetEntriesFast()>0){
 511     prevPHOS->AddFirst(fPHOSEvent) ;
 512     fPHOSEvent=0;
 513     if(prevPHOS->GetSize()>100){//Remove redundant events
 514       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
 515       prevPHOS->RemoveLast() ;
 516       delete tmp ;
 517     }
 518   }
 519
 520   // Post output data.
 521   PostData(1, fOutputContainer);
 522   fEventCounter++;
 523 }
 524
 525 //________________________________________________________________________
 526 void AliAnalysisTaskPi0Calib::Terminate(Option_t *)
 527 {
 528   // Draw result to the screen
 529   // Called once at the end of the query
 530
 531 }
 532
 533 //________________________________________________________________________
 534 Bool_t AliAnalysisTaskPi0Calib::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
 535 {
 536   //Check if this channel belogs to the good ones
 537
 538   if(strcmp(det,"PHOS")==0){
 539     if(mod>5 || mod<1){
 540       AliError(Form("No bad map for PHOS module %d ",mod)) ;
 541       return kTRUE ;
 542     }
 543     if(!fPHOSBadMap[mod]){
 544       AliError(Form("No Bad map for PHOS module %d",mod)) ;
 545       return kTRUE ;
 546     }
 547     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
 548       return kFALSE ;
 549     else
 550       return kTRUE ;
 551   }
 552   else{
 553     AliError(Form("Can not find bad channels for detector %s ",det)) ;
 554   }
 555   return kTRUE ;
 556 }
 557 //_____________________________________________________________________________
 558 void AliAnalysisTaskPi0Calib::FillHistogram(const char * key,Double_t x)const{
 559   //FillHistogram
 560   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
 561   if(tmpI){
 562     tmpI->Fill(x) ;
 563     return ;
 564   }
 565   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
 566   if(tmpF){
 567     tmpF->Fill(x) ;
 568     return ;
 569   }
 570   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
 571   if(tmpD){
 572     tmpD->Fill(x) ;
 573     return ;
 574   }
 575   AliInfo(Form("can not find histogram <%s> ",key)) ;
 576 }
 577 //_____________________________________________________________________________
 578 void AliAnalysisTaskPi0Calib::FillHistogram(const char * key,Double_t x,Double_t y)const{
 579   //FillHistogram
 580   TObject * tmp = fOutputContainer->FindObject(key) ;
 581   if(!tmp){
 582     AliInfo(Form("can not find histogram <%s> ",key)) ;
 583     return ;
 584   }
 585   if(tmp->IsA() == TClass::GetClass("TH1F")){
 586     ((TH1F*)tmp)->Fill(x,y) ;
 587     return ;
 588   }
 589   if(tmp->IsA() == TClass::GetClass("TH2F")){
 590     ((TH2F*)tmp)->Fill(x,y) ;
 591     return ;
 592   }
 593   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
 594 }
 595
 596 //_____________________________________________________________________________
 597 void AliAnalysisTaskPi0Calib::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
 598   //Fills 1D histograms with key
 599   TObject * tmp = fOutputContainer->FindObject(key) ;
 600   if(!tmp){
 601     AliInfo(Form("can not find histogram <%s> ",key)) ;
 602     return ;
 603   }
 604   if(tmp->IsA() == TClass::GetClass("TH2F")){
 605     ((TH2F*)tmp)->Fill(x,y,z) ;
 606     return ;
 607   }
 608   if(tmp->IsA() == TClass::GetClass("TH3F")){
 609     ((TH3F*)tmp)->Fill(x,y,z) ;
 610     return ;
 611   }
 612 }
 613 //_____________________________________________________________________________
 614 Bool_t AliAnalysisTaskPi0Calib::TestLambda(Double_t pt,Double_t l1,Double_t l2){
 615
 616   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 617   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 618   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 619   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 620   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 621   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 622               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 623               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 624   return (R2<2.5*2.5) ;
 625
 626 }
 627 //_____________________________________________________________________________
 628 Bool_t AliAnalysisTaskPi0Calib::TestTOF(Double_t /*e*/, Double_t tof){
 629
 630   Double_t res = 100.e-9 ; //3.31714e-09  + 1.42519e-08*TMath::Exp(-e/4.61728e-01) ;
 631   Double_t mean=0. ;
 632
 633   return TMath::Abs(tof-mean)<res ;
 634
 635 }
 636 //____________________________________________________________________________
 637 Double_t AliAnalysisTaskPi0Calib::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
 638   //Parameterization of LHC10h period
 639   //_true if neutral_
 640
 641   Double_t meanX=0;
 642   Double_t meanZ=0.;
 643   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
 644               6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
 645   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
 646   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 647   Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
 648
 649   if(mf<0.){ //field --
 650     meanZ = -0.468318 ;
 651     if(charge>0)
 652       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
 653     else
 654       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
 655   }
 656   else{ //Field ++
 657     meanZ= -0.468318;
 658     if(charge>0)
 659       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
 660     else
 661       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
 662   }
 663
 664   Double_t rz=(dz-meanZ)/sz ;
 665   Double_t rx=(dx-meanX)/sx ;
 666   return TMath::Sqrt(rx*rx+rz*rz) ;
 667 }
 668 //____________________________________________________________________________