PWG1/ITS/AliAnalysisTaskdEdxSSDQA.cxx

   1 //SSD dEdx gain calibration task
   2 //to run in QA train
   3 //autor Marek Chojnacki
   4 //Marek.Chojnacki@cern.ch
   5 // last line
   6 #include "TChain.h"
   7 #include "TTree.h"
   8 #include "TH1F.h"
   9 #include "TH2F.h"
  10 #include "TH3F.h"
  11 #include "TFile.h"
  12
  13 #include "AliAnalysisTaskdEdxSSDQA.h"
  14 #include "AliAnalysisManager.h"
  15
  16 #include "AliESDEvent.h"
  17 #include "AliESDInputHandler.h"
  18
  19
  20 #include "AliESDtrack.h"
  21 #include "AliExternalTrackParam.h"
  22 #include "Riostream.h"
  23
  24 #include "AliESDfriend.h"
  25 #include "AliESDfriendTrack.h"
  26 #include "AliTrackPointArray.h"
  27 #include "AliCDBManager.h"
  28 #include "AliGeomManager.h"
  29 #include "TMath.h"
  30 #include "TGeoMatrix.h"
  31
  32 using namespace std;
  33 ClassImp(AliAnalysisTaskdEdxSSDQA)
  34
  35 //________________________________________________________________________
  36 AliAnalysisTaskdEdxSSDQA::AliAnalysisTaskdEdxSSDQA(const char *name)
  37 : AliAnalysisTaskSE(name), fHist1(0),fHist2(0),fHist3(0),fHist4(0),fHist5(0),fHist6(0),fListOfHistos(0),fPcut(0.0),fdothecorrection(0)
  38 {
  39   // Constructor
  40         // Define input and output slots here
  41   // Input slot #1 works with a TChain
  42  // fcorrections=new Float_t[1698*12];
  43   for(int i=0;i<1698*12;i++)
  44         fcorrections[i]=1.0;
  45   DefineOutput(1, TList::Class());
  46 }
  47
  48 //________________________________________________________________________
  49 void AliAnalysisTaskdEdxSSDQA::UserCreateOutputObjects()
  50 {
  51   // Create histograms
  52   // Called once
  53
  54         fListOfHistos = new TList();
  55         fHist1 =new TH2F("QAChargeRatio","QAChargeRatio;Module;CR",1698,-0.5,1697.5,80,-1.0,1.0);
  56         fListOfHistos->Add(fHist1);
  57         fHist2=new TH2F("QACharge","QACharge;Module;Q",1698,-0.5,1697.5,150,0,300);
  58         fListOfHistos->Add(fHist2);
  59         fHist3=new TH3F("ChargeRatiovCharge","CRvQ;module;Q;CR",1698,-0.5,1697.5,80,0,1600,50,-1.0,1.0);
  60         fListOfHistos->Add(fHist3);
  61         fHist4=new TH2F("QAChargedinchips","chargedinchips;chip;Q",1698*12,-0.5,1698*12-0.5,150,0,300);
  62         fListOfHistos->Add(fHist4);
  63         if(fdothecorrection)
  64         {
  65                 fHist5=new TH2F("QAChargeCorrected","chargedCorrected;Module;QCorrected",1698,-0.5,1697.5,150,0,300);
  66                 fListOfHistos->Add(fHist5);
  67         }
  68         fHist6=new TH2F("QNvQP",";QN;QP",160,0,1600,1600,0,1600);
  69         fListOfHistos->Add(fHist6);
  70
  71         PostData(1,  fListOfHistos);
  72 }
  73 //________________________________________________________________________
  74 void AliAnalysisTaskdEdxSSDQA::LocalInit()
  75 {
  76
  77         Printf("end of LocalInit");
  78 }
  79 //______________________________________________________________________
  80 AliAnalysisTaskdEdxSSDQA::~AliAnalysisTaskdEdxSSDQA()
  81 {
  82 }
  83
  84 //________________________________________________________________________
  85
  86 void AliAnalysisTaskdEdxSSDQA::UserExec(Option_t *)
  87 {
  88     // Main loop
  89     // Called for each event
  90     AliESDEvent* esd = dynamic_cast<AliESDEvent*> (InputEvent());
  91     if(!esd)
  92     {
  93         Printf("ERROR: Input ESD Event not available");
  94         PostData(1,  fListOfHistos);
  95         return;
  96     }
  97
  98     Int_t nTracks=esd->GetNumberOfTracks();
  99
 100
 101
 102         if(!ESDfriend())
 103         {
 104                 Printf("problem with friend");
 105                 PostData(1,  fListOfHistos);
 106                 return;
 107         }
 108         //Printf("Event nTracks %d", nTracks);
 109         AliTrackPointArray*  trackar=0x0;
 110         Bool_t l5;
 111         Bool_t l6;
 112         Int_t npoints;
 113         AliTrackPoint point;
 114         Int_t nlayer=0;
 115         Int_t id=0;
 116         Float_t chargeratio=0.0;
 117         Float_t charge=0.0;
 118         for(int itr=0;itr<nTracks;itr++)
 119         {
 120                 AliESDtrack* track= esd->GetTrack(itr);
 121
 122                 if(TMath::Abs(track->Eta())>0.9)
 123                         continue;
 124                 if(track->GetP()>10.0)
 125                         continue;
 126                 if (track->IsOn(AliESDtrack::kITSrefit)&&track->IsOn(AliESDtrack::kTPCrefit))
 127                 {
 128                          l5=track->HasPointOnITSLayer(4);
 129                          l6=track->HasPointOnITSLayer(5);
 130                         if (!(l5||l6))//only tracks with SSD point
 131                                 continue;
 132                         Double_t tmpQESD[4]={-1.0,-1.0,-1.0,-1.0};
 133                         track->GetITSdEdxSamples(tmpQESD);
 134                         trackar=(AliTrackPointArray*)track->GetTrackPointArray();
 135                         if(!trackar)
 136                                 continue;
 137                         npoints=trackar->GetNPoints();
 138                         for(int itnp=0;itnp<npoints;itnp++)
 139                         {
 140                                 if(trackar->GetPoint(point,itnp))
 141                                 {
 142                                         nlayer=AliGeomManager::VolUIDToLayer(point.GetVolumeID(), id);//layer number
 143                                 }
 144                                 else
 145                                         continue;
 146                                 if(nlayer==5||nlayer==6)
 147                                 {
 148                                         TGeoHMatrix* geomatrix=AliGeomManager::GetMatrix(point.GetVolumeID());
 149                                         //geomatrix->Print();
 150                                         if(point.GetCharge()>0.0&&point.IsExtra()==kFALSE)//do not use additional clusters
 151                                         {
 152                                                 Double_t local[3]={0.0,0.0,0.0};
 153                                                 Double_t global[3]={0.0,0.0,0.0};
 154                                                 global[0]=point.GetX();
 155                                                 global[1]=point.GetY();
 156                                                 global[2]=point.GetZ();
 157                                                 geomatrix->MasterToLocal(global,local);
 158                                                  chargeratio=point.GetChargeRatio();
 159                                                  charge=point.GetCharge();
 160                                                  Float_t fQNnotcorr=charge*(1.0+chargeratio);
 161                                                  Float_t fQPnotcorr=charge*(1.0-chargeratio);
 162                                                  fHist6->Fill(fQNnotcorr,fQPnotcorr);
 163                                                 // cout<<point.GetCharge()<<" "<<point.GetChargeRatio()<<" "<<nlayer<<" "<<id<<endl;
 164                                                 if(nlayer==5&&tmpQESD[2]>0.0)
 165                                                 {
 166                                                         fHist1->Fill(id,chargeratio);
 167                                                         if(track->GetP()>fPcut)
 168                                                         {
 169                                                                 Float_t fQP=tmpQESD[2]*(1.0-chargeratio);
 170                                                                 Float_t fQN=tmpQESD[2]*(1.0+chargeratio);
 171                                                                 Int_t fPchip=Pstrip5(10000.0*local[0],10000.0*local[2])/128;
 172                                                                 Int_t fNchip=Nstrip5(10000.0*local[0],10000.0*local[2])/128;
 173
 174                                                                 fHist2->Fill(id,tmpQESD[2]);
 175                                                                 fHist4->Fill(id*12+fPchip,fQP);
 176                                                                 fHist4->Fill(id*12+6+fNchip,fQN);
 177                                                                 if(fdothecorrection)
 178                                                                         fHist5->Fill(id,0.5*(fQP*fcorrections[12*id+fPchip]+fQN*fcorrections[12*id+6+fNchip]));
 179                                                         }
 180                                                         fHist3->Fill(id,charge,chargeratio);
 181                                                 }
 182                                                 if(nlayer==6&&tmpQESD[3]>0.0)
 183                                                 {
 184                                                         fHist1->Fill(id+748,chargeratio);
 185                                                         if(track->GetP()>fPcut)
 186                                                         {
 187
 188                                                                 Float_t fQP=tmpQESD[3]*(1.0-chargeratio);
 189                                                                 Float_t fQN=tmpQESD[3]*(1.0+chargeratio);
 190                                                                 Int_t fPchip=Pstrip6(10000.0*local[0],10000.0*local[2])/128;
 191                                                                 Int_t fNchip=Nstrip6(10000.0*local[0],10000.0*local[2])/128;
 192
 193                                                                 fHist2->Fill(id+748,tmpQESD[3]);
 194                                                                 fHist4->Fill((id+748)*12+fPchip,fQP);
 195                                                                 fHist4->Fill((id+748)*12+6+fNchip,fQN);
 196                                                                 if(fdothecorrection)
 197                                                                         fHist5->Fill((id+748),0.5*(fQP*fcorrections[12*(id+748)+fPchip]+fQN*fcorrections[12*(id+748)+6+fNchip]));
 198
 199                                                         }
 200                                                         fHist3->Fill(id+748,charge,chargeratio);
 201                                                 }
 202
 203                                         }
 204                                 }
 205                         }
 206                 }
 207         }
 208         //Printf("Event nTracks end %d", nTracks);
 209         // Post output data.
 210         PostData(1,  fListOfHistos);
 211 }
 212
 213 //________________________________________________________________________
 214 void AliAnalysisTaskdEdxSSDQA::Terminate(Option_t *)
 215 {
 216         //terminate
 217         Printf("end of Terminate");
 218 }
 219 //_____________________________________________________________________________
 220 Int_t  AliAnalysisTaskdEdxSSDQA::Pstrip5(Float_t x,Float_t z) const
 221 {
 222         // P strip for layer 5 from local cooridnates
 223         Float_t value;
 224         x=x+36432.5-223.9;
 225         z=z-(76000.0/7.0);
 226         value=-(3.0/38000.0)*z+x/95.0;
 227         return TMath::Nint(value);
 228 }
 229 //___________________________________________________________________________
 230 Int_t  AliAnalysisTaskdEdxSSDQA::Nstrip5(Float_t x,Float_t z)  const
 231 {
 232         // N strip for layer 5 from local cooridnates
 233         Float_t value;
 234         x=x+36432.5-223.9;
 235         z=z-(76000.0/7.0);
 236         value=x/95.0+(11.0/38000.0)*z;
 237         return TMath::Nint(value);
 238 }
 239 //___________________________________________________________________________
 240 Int_t  AliAnalysisTaskdEdxSSDQA::Pstrip6(Float_t x,Float_t z)  const
 241 {
 242         // P strip for layer 6 from local cooridnates
 243         Float_t value;
 244         x=x-36432.5-223.9;
 245         z=z+(76000.0/7.0);
 246         value=-(3.0/38000.0)*z-x/95.0;
 247         return TMath::Nint(value);
 248 }
 249 //__________________________________________________________________________________
 250 Int_t  AliAnalysisTaskdEdxSSDQA::Nstrip6(Float_t x,Float_t z)  const
 251 {
 252         // N strip for layer 6 from local cooridnates
 253         Float_t value;
 254         x=x-36432.5-223.9;
 255         z=z+(76000.0/7.0);
 256         value=x*(-1.0/95.0)+(11.0/38000.0)*z;
 257         return TMath::Nint(value);
 258 }
 259 //___________________________________________________________________________________
 260 void  AliAnalysisTaskdEdxSSDQA::SetDoChipCorretions(const char* filename)
 261 {
 262 //Upload corrections for each chip only for test
 263         cout<<filename<<endl;
 264         ifstream infile(filename);
 265         if(!infile.good())
 266                 return;
 267         Float_t value=0.0;
 268         fdothecorrection=1;
 269         for (int i=0;i<1698;i++)
 270         {
 271                 for (int j=0;j<6;j++)
 272                 {
 273                         infile>>value;
 274                         cout<<value<<" ";
 275                         if(value>0.0)
 276                                 fcorrections[i*12+j]=84.0/value;
 277                         else
 278                                 fcorrections[i*12+j]=1.0;
 279                         infile>>value;
 280                         infile>>value;
 281                         cout<<value<<" ";
 282                         if(value>0.0)
 283                                 fcorrections[i*12+j+6]=84.0/value;
 284                         else
 285                                 fcorrections[i*12+j+6]=1.0;
 286                         infile>>value;
 287                 }
 288                 cout<<endl;
 289         }
 290 }
 291