ITS/ITSSPDVertexDiamondda.cxx

   1 /*
   2 Contact: cvetan.cheshkov@cern.ch
   3 Link: missing
   4 Run Type: PHYSICS
   5 DA Type: MON
   6 Number of events needed: 10000
   7 Input Files:
   8 Output Files:
   9 Trigger types used: PHYSICS
  10 */
  11
  12 #define X_LIMIT 2.0
  13 #define Y_LIMIT 2.0
  14 #define Z_LIMIT 50.0
  15 #define X_DELTA 0.02
  16 #define Y_DELTA 0.02
  17 #define Z_DELTA 0.2
  18 #define OUTPUT_FILE "SPDVertexDiamondDA.root"
  19 #define CDB_STORAGE "local://$ALICE_ROOT"
  20 #define N_EVENTS_AUTOSAVE 50
  21 #define NFITPARAMS 5
  22
  23 extern "C" {
  24 #include "daqDA.h"
  25 }
  26
  27 #include "event.h"
  28 #include "monitor.h"
  29
  30 #include <TTree.h>
  31 #include <TH1.h>
  32 #include <TH2.h>
  33 #include <TF2.h>
  34 #include <TFile.h>
  35 #include <TPluginManager.h>
  36 #include <TROOT.h>
  37 #include <TFitter.h>
  38
  39 #include "AliRawReaderDate.h"
  40 #include "AliGeomManager.h"
  41 #include "AliCDBManager.h"
  42 #include "AliESDVertex.h"
  43 #include "AliITSDetTypeRec.h"
  44 #include "AliITSInitGeometry.h"
  45 #include "AliITSVertexer3DTapan.h"
  46
  47 TF2 *fitFcn = 0x0;
  48
  49 AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ);
  50
  51 int main(int argc, char **argv) {
  52
  53   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
  54                                         "*",
  55                                         "TStreamerInfo",
  56                                         "RIO",
  57                                         "TStreamerInfo()");
  58   TFitter *minuitFit = new TFitter(NFITPARAMS);
  59   TVirtualFitter::SetFitter(minuitFit);
  60
  61   int status;
  62   if (argc<2) {
  63     printf("Wrong number of arguments\n");
  64     return -1;
  65   }
  66
  67   /* define data source : this is argument 1 */
  68   status=monitorSetDataSource( argv[1] );
  69   if (status!=0) {
  70     printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
  71     return -1;
  72   }
  73
  74   /* declare monitoring program */
  75   status=monitorDeclareMp( __FILE__ );
  76   if (status!=0) {
  77     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
  78     return -1;
  79   }
  80
  81   /* define wait event timeout - 1s max */
  82   monitorSetNowait();
  83   monitorSetNoWaitNetworkTimeout(1000);
  84
  85   /* log start of process */
  86   printf("Vertex-Diamond SPD DA started\n");
  87
  88   /* init some counters */
  89   int nevents_with_vertex = 0;
  90   int nevents_physics=0;
  91   int nevents_total=0;
  92
  93   struct eventHeaderStruct *event;
  94   eventTypeType eventT;
  95
  96   // Get run number
  97   if (getenv("DATE_RUN_NUMBER")==0) {
  98     printf("DATE_RUN_NUMBER not properly set.\n");
  99     return -1;
 100   }
 101   int runNr = atoi(getenv("DATE_RUN_NUMBER"));
 102
 103   // Histograms initialization
 104   TH2F *hXY = new TH2F("hXY","Vertex Diamond (Y vs X)",
 105                        2*(Int_t)(X_LIMIT/X_DELTA),-X_LIMIT,X_LIMIT,
 106                        2*(Int_t)(Y_LIMIT/Y_DELTA),-Y_LIMIT,Y_LIMIT);
 107   TH1F *hZ  = new TH1F("hZ"," Longitudinal Vertex Profile",
 108                        2*(Int_t)(Z_LIMIT/Z_DELTA),-Z_LIMIT,Z_LIMIT);
 109
 110   // Global initializations
 111   AliCDBManager *man = AliCDBManager::Instance();
 112   man->SetDefaultStorage(CDB_STORAGE);
 113   man->SetRun(runNr);
 114   AliGeomManager::LoadGeometry("geometry.root");
 115   AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 116
 117   // ITS initializations
 118   AliITSInitGeometry initgeom;
 119   AliITSgeom *geom = initgeom.CreateAliITSgeom();
 120   printf("Geometry name: %s\n",(initgeom.GetGeometryName()).Data());
 121
 122   AliITSDetTypeRec *detTypeRec = new AliITSDetTypeRec();
 123   detTypeRec->SetITSgeom(geom);
 124   detTypeRec->SetDefaults();
 125   detTypeRec->SetDefaultClusterFindersV2(kTRUE);
 126
 127   /* main loop (infinite) */
 128   for(;;) {
 129
 130     /* check shutdown condition */
 131     if (daqDA_checkShutdown()) {break;}
 132
 133     /* get next event (blocking call until timeout) */
 134     status=monitorGetEventDynamic((void **)&event);
 135     if (status==MON_ERR_EOF) {
 136       printf ("End of File detected\n");
 137       break; /* end of monitoring file has been reached */
 138     }
 139
 140     if (status!=0) {
 141       printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 142       break;
 143     }
 144
 145     /* retry if got no event */
 146     if (event==NULL) {
 147       continue;
 148     }
 149
 150     nevents_total++;
 151     eventT=event->eventType;
 152     switch (event->eventType){
 153
 154       /* START OF RUN */
 155     case START_OF_RUN:
 156       break;
 157       /* END START OF RUN */
 158
 159     /* END OF RUN */
 160     case END_OF_RUN:
 161       break;
 162
 163     case PHYSICS_EVENT:
 164       nevents_physics++;
 165       AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 166
 167       // Run SPD cluster finder
 168       TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
 169       detTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
 170
 171       // Run vertex-finder
 172       AliITSVertexer3DTapan *vertexer = new AliITSVertexer3DTapan(geom,1000);
 173       vertexer->LoadClusters(clustersTree);
 174       AliESDVertex *vtx = new AliESDVertex();
 175       vertexer->FindVertexForCurrentEvent(vtx);
 176
 177       if (TMath::Abs(vtx->GetChi2()) < 0.1) {
 178         // Fill the vertex into the histos
 179         nevents_with_vertex++;
 180         hXY->Fill(vtx->GetXv(),vtx->GetYv());
 181         hZ->Fill(vtx->GetZv());
 182
 183         // Auto save
 184         if ((nevents_with_vertex%N_EVENTS_AUTOSAVE) == 0) {
 185           TFile outFile(OUTPUT_FILE, "update");
 186           AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
 187           if (fitVtx) {
 188             fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
 189             TH1 *fithXY = fitFcn->CreateHistogram();
 190             fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
 191             delete fithXY;
 192           }
 193           hXY->Write(hXY->GetName(),TObject::kOverwrite);
 194           hZ->Write(hZ->GetName(),TObject::kOverwrite);
 195           outFile.Close();
 196           delete fitVtx;
 197         }
 198       }
 199
 200       delete vtx;
 201       delete vertexer;
 202
 203       delete clustersTree;
 204       delete rawReader;
 205
 206     }
 207
 208     /* free resources */
 209     free(event);
 210
 211     /* exit when last event received, no need to wait for TERM signal */
 212     if (eventT==END_OF_RUN) {
 213       printf("EOR event detected\n");
 214       break;
 215     }
 216   }
 217
 218   if (detTypeRec) delete detTypeRec;
 219
 220   // Store the final histograms
 221   TFile outFile(OUTPUT_FILE, "update");
 222   if (nevents_with_vertex > N_EVENTS_AUTOSAVE) {
 223     // Fit XY & Z histograms
 224     AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
 225     if (fitVtx) {
 226       fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
 227       TH1 *fithXY = fitFcn->CreateHistogram();
 228       fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
 229       delete fithXY;
 230     }
 231     delete fitVtx;
 232   }
 233   hXY->Write(hXY->GetName(),TObject::kOverwrite);
 234   hZ->Write(hZ->GetName(),TObject::kOverwrite);
 235   outFile.Close();
 236
 237   delete minuitFit;
 238   TVirtualFitter::SetFitter(0);
 239
 240   /* write report */
 241   printf("Run #%s, received %d events with vertex, out of %d physics and out of %d total events\n",getenv("DATE_RUN_NUMBER"),nevents_with_vertex,nevents_physics,nevents_total);
 242
 243   status=0;
 244
 245   /* export file to FXS */
 246   if (daqDA_FES_storeFile(OUTPUT_FILE, "VertexDiamond")) {
 247     status=-2;
 248   }
 249
 250   return status;
 251 }
 252
 253 Double_t fitFunction(Double_t *x, Double_t *par) {
 254
 255   Double_t t1 =   x[0] - par[1];
 256   Double_t t2 =   x[1] - par[2];
 257
 258   return par[0]*TMath::Exp(-0.5*(t1*t1/(par[3]*par[3])+t2*t2/(par[4]*par[4])));
 259 }
 260
 261 AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ)
 262 {
 263
 264   if (!fitFcn) {
 265     fitFcn = new TF2("fitFcn",fitFunction,
 266                      -X_LIMIT,X_LIMIT,
 267                      -Y_LIMIT,Y_LIMIT,NFITPARAMS);
 268     fitFcn->SetNpx(2*(Int_t)(X_LIMIT/X_DELTA));
 269     fitFcn->SetNpy(2*(Int_t)(Y_LIMIT/Y_DELTA));
 270     fitFcn->SetParameters(hXY->GetMaximum(),0,0,hXY->GetRMS(1),hXY->GetRMS(2));
 271     fitFcn->Update();
 272   }
 273
 274   if (hXY->Fit("fitFcn","L0V+") != 0) {
 275     printf("XY fit failed!");
 276     return 0x0;
 277   }
 278
 279   Double_t pos[3],poserr[3];
 280   pos[0] = fitFcn->GetParameter(1);
 281   pos[1] = fitFcn->GetParameter(2);
 282   poserr[0] = fitFcn->GetParameter(3);
 283   poserr[1] = fitFcn->GetParameter(4);
 284
 285   // Could be replaced by something more robust...
 286   pos[2] = hZ->GetMean();
 287   poserr[2] = hZ->GetRMS();
 288
 289   return new AliESDVertex(pos,poserr);
 290 }